# AGN energetics and lifetimes from remnant radio galaxies

Benjamin Quici1, Ross J. Turner2, Nicholas Seymour1, Natasha Hurley-Walker1
1International Centre for Radio Astronomy Research, Curtin University,
Bentley, WA 6102, Australia
2School of Natural Sciences, University of Tasmania, Private Bag 37, Hobart,
7001, Australia<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The energy coupling efficiency of active galactic nucleus (AGN) outbursts is
known to differ significantly with factors including the jet kinetic power,
duration of the outburst, and properties of the host galaxy cluster. As such,
constraints on their jet power and lifetime functions are crucial to quantify
the role of kinetic-mode AGN feedback on the evolution of galaxies since
$z\sim 1$. In this work, we address this issue by measuring the energetics of
a sample of 79 low-redshift (0.02 $<z<$ 0.2) remnant radio galaxies compiled
from large-sky radio surveys – these objects uniquely probe the full duration
of an AGN outburst. The jet kinetic power and outburst duration of each
remnant are determined using the RAiSE dynamical model based on the surface
brightness distribution observed in multi-frequency radio images. We compare
the energetics constrained for this sample to those predicted for mock radio
source populations – with various intrinsic functions for jet power and
lifetime distributions – to correct for sample selection biases imposed on our
sample. The intrinsic jet power and lifetime functions that yield a selection-
biased mock population most similar to our observed sample are found using
Bayesian inference. Our analysis places robust constraints on assumed power-
law indices for the intrinsic jet power and lifetime functions: $p(Q)\propto
Q^{-1.49\pm 0.07}$ and $p(t_{\mathrm{on}})\propto t_{\mathrm{on}}^{-0.97\pm
0.12}$ respectively. We discuss the implication of these findings for
feedback-regulated accretion and the self-regulating nature of jet activity.
The methodology proposed in this work can be extended to active radio galaxies
in future studies.

###### keywords:

galaxies: active – galaxies: jets – radio continuum: galaxies.

††pubyear: 2024††pagerange: AGN energetics and lifetimes from remnant radio
galaxies–AGN energetics and lifetimes from remnant radio galaxies

##  1 Introduction

In the concordance view of galaxy formation and evolution, relativistic jets
associated with radio-loud active galactic nuclei (AGNs) play a critical role
in regulating the gas residing in the interstellar and circumgalactic
environments. On small (galactic) scales, the jets are capable of dispelling
large amounts of molecular gas (e.g., Nesvadba et al., 2008; Dasyra & Combes,
2011; Morganti et al., 2013; Emonts et al., 2014; Morganti et al., 2015;
Kakkad et al., 2018; Morganti et al., 2021), and can seemingly promote
(Nesvadba et al., 2008; Dasyra & Combes, 2012; Morganti et al., 2013;
Mukherjee et al., 2021) or suppress (Croft et al., 2006; Tortora et al., 2009;
Crockett et al., 2012; Dugan et al., 2017) star formation. On larger scales,
the buoyant rise of jet-inflated bubbles excavate cavities in their hot
intracluster environments (e.g., Churazov et al., 2002), which act to suppress
the cooling flows at the centre of cool-core clusters (Fabian et al., 2003;
Forman et al., 2005; McNamara & Nulsen, 2007; Mittal et al., 2009; Alexander &
Hickox, 2012; Fabian, 2012). Recent findings also show that ‘maintenance-mode
feedback’, the self-regulating mechanism by which the jets maintain a
heating/cooling balance in their hot atmospheres, must be invoked in galaxy
formation models to explain, in particular, the suppression of star formation
in the most massive galaxies (e.g., Bower et al., 2006; Croton et al., 2006;
Shabala & Alexander, 2009; Raouf et al., 2017).

Despite the various successes of modern galaxy formation models (e.g., Dubois
et al., 2014; Vogelsberger et al., 2014; Schaye et al., 2015), constraints
from observations are still required to quantify the global energetic impact
of AGNs, and to elucidate the connection between the environment and fuelling
of the supermassive black hole (SMBH). One valuable insight lies in the
observed mass-dependence of the radio loud fraction, $f_{\mathrm{RL}}$, which
rises with the stellar mass of the host galaxy (Best et al., 2005) and likely
achieves ubiquity ($\approx 100\%$) amongst the most massive galaxies (Sabater
et al., 2019); precisely the objects in which feedback is required the most.
Treating $f_{\mathrm{RL}}$ as a proxy for the AGN jet duty cycle, the result
implies a clear connection between the fuelling of the SMBH and the
surrounding environment. Yet, questions still remain regarding the nature of
this relationship. Principally, if it is AGN feedback that links the quenching
of star formation to the prevalence of the jets, then the energies involved
must be quantified to understand how this process occurs. Additionally, the
$f_{\mathrm{RL}}$ cannot separate between objects with rapid versus long-
period duty-cycles; this distinction is important, considering that such
objects will have vastly different duration active jet phases, and will thus
differ in their feedback efficiencies (Yates et al., 2018).

Valuable constraints towards AGN feedback mechanisms can be derived by
examining the distribution of AGN jet kinetic powers and lifetimes. For a
single AGN outburst, its total energy output can be found by integrating the
instantaneous jet kinetic power (the jet power herein) over the the duration
of the outburst (the jet lifetime herein). Considering that the jets are
ultimately produced by accretion onto the SMBH (Lynden-Bell, 1969), we can
expect each of these jet properties to be set by the conditions of and around
the black hole; namely mass accretion rate and spin for jet power (Daly, 2009,
2016, 2021), and the underlying accretion mechanism for the jet outburst
(Novak et al., 2011; Gabor & Bournaud, 2013). It follows that the distribution
in jet powers and lifetimes not only constrains the energies released by AGNs,
but also begins to probe the underlying mechanisms responsible for the
production and fuelling of the jets. Considering also the mass-dependence of
the AGN jet duty cycle, reconciling the global properties of the jets together
with those of their environment are an important step towards linking AGN
fuelling and feedback.

Throughout the decades, a variety of techniques have been presented to measure
the jet powers and lifetimes of radio-loud AGNs. The diffuse lobes inflated by
the jets contain synchrotron-emitting plasma that forms an approximately
power-law spectrum at radio frequencies. The energy-dependent synchrotron loss
rate allows for a derivation of the plasma age based on the observed spectral
steepening in the radio spectrum (e.g., Parma et al., 2007; Murgia et al.,
2011; Brienza et al., 2016; Duchesne et al., 2020), and is referred to in the
literature as the ‘spectral ageing’ method (e.g., Myers & Spangler, 1985). A
shortcoming of this method is the requirement of the lobe magnetic field
strength, which in turn requires knowledge of the equipartition factor and the
pressure in the ambient environment. By making several assumptions regarding
the lobes, e.g., their composition and minimum energy, inferences of the jet
power can be made based on their observed radio luminosity (Willott et al.,
1999). Comparisons of these estimates against more sophisticated methods shows
a general agreement (e.g., see Ineson et al., 2017; Turner et al., 2018b),
however with an appreciable level of scatter in part due to the departure from
the underlying assumptions, but also due to the confounding impact of the
source age and environment (Godfrey & Shabala, 2013; Hardcastle & Krause,
2013). Cavities in the hot X-ray emitting intracluster gas allow for robust
estimates to be made of the $pdV$ work done by the expanding lobes; with an
assumption about the source age (discussed above), inferences of the jet power
can then be made. However, this technique is generally limited to low
redshifts ($z\lesssim 0.5$) as the trade-off between resolution and surface
brightness sensitivity in X-ray observations prevents measurements of the
shape of the gas density profile beyond the cluster core (e.g., Turner &
Shabala, 2020).

Dynamical model-based methods have increasingly proved an attractive
alternative for measuring the intrinsic properties of AGN jets. The current
generation of analytical models (e.g.,Turner et al. 2018a; Hardcastle 2018;
Turner et al. 2023; see Turner & Shabala 2023 for a review) match the
evolutionary histories predicted by modern relativistic hydrodynamic
simulations (e.g., Hardcastle & Krause, 2013, 2014; English et al., 2016;
Yates et al., 2018; English et al., 2019; Yates-Jones et al., 2021) whilst
including a full treatment of radiative and expansion losses in their
synchrotron emissivity calculations. The energetics of radio-loud AGNs can be
constrained by comparing simulated attributes of synthetic radio sources
(e.g., size and luminosity) against those for the observed objects. This
comparison is typically done through a Bayesian parameter inversion (e.g., see
Turner & Shabala, 2015; Turner et al., 2018b; Turner, 2018; Mahatma et al.,
2020; Turner et al., 2020; Quici et al., 2022), allowing the intrinsic radio
source parameters to be fitted. Such methods have the ability to
simultaneously constrain the jet power and active lifetime of a radio source,
in addition to the lobe magnetic field strength (Turner et al., 2018b) or host
redshift (Turner et al., 2020), based on its angular size and multi-frequency
radio spectrum. Quici et al. (2022) showed that high-resolution radio images
at two frequencies (e.g., 151 MHz and 1.4 GHz; typical of large-sky surveys)
can provide sufficient constraining power in a Bayesian parameter inversion to
fit at least four intrinsic jet parameters including the duration of the
remnant phase, $t_{\rm rem}$. However, such parameter inversions rely on the
ability of the underlying model to capture the necessary physics needed to
explain the observation. For example, the analytical model of Turner et al.
(2023) does not currently capture processes such as jet precession or a
bending of the jet, and may therefore be inapplicable for certain active radio
sources; on the other hand, these issues are less critical for the remnant
phase as regions of enhanced emission associated with the jets will rapidly
fade due to adiabatic expansion and radiative loss mechanisms (e.g., see
Mahatma 2023 for a review).

Recently, Jurlin et al. (2020) selected an observationally-complete sample of
radio galaxies from the LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al.,
2019, 2022). Through the application of forward modelling facilitated by the
Radio AGN in Semi-analytic Environments (RAiSE; see Turner et al., 2023)
model, Shabala et al. (2020) used this sample to constrain the intrinsic (or
true) jet power and lifetime functions (see also Hardcastle et al. 2019 for a
similar approach). Their approach (i.e., that of Shabala et al., 2020) relied
on the comparison between mock and observed distributions in flux density and
angular size, as well as the constraints from the observed fraction of remnant
and restarted radio sources. Their work demonstrated that power-law lifetime
models were needed to reproduce the observed properties of their reference
sample, unlike constant-age models which under-predict the fraction of remnant
radio galaxies.

In this work, we combine and extend the approaches taken by Shabala et al.
(2020) and Quici et al. (2022) to quantify the intrinsic jet power and
lifetime functions for a sample of remnant radio galaxies. Such objects are
unique in that they allow for a direct measurement of the full duration of the
jet outburst, thus offering tight constraints on the AGN jet lifetime
function. First, we identify a sample of 79 remnant radio galaxies brighter
than $0.5$ Jy, with angular sizes $\theta\geqslant 4^{\prime}$, and within a
redshift range of $0.02<z<0.2$ (section 2). Next, we adopt the method of Quici
et al. (2022) to fit their intrinsic parameters, allowing us to probe their
observed jet power and active age distributions (section 3). To deal with the
inherent selection biases, we construct mock remnant catalogues based on an
assumed jet power and lifetime function, and filter the resulting radio source
population by our sample selection criteria (section 4). We test a grid of jet
power and lifetime models, and compare their predictions to our observed
sample statistics to constrain the true jet power and lifetime functions
(section 5). We discuss our results in the context of the broader literature
and conclude in section 6.

Throughout this work, we assume the $\Lambda$CDM concordance cosmology with
$\Omega_{M}$ = 0.3, $\Omega_{\Lambda}=0.7$ and $H_{0}$ = 70 km s-1 Mpc-1
(Planck Collaboration et al., 2016), and the J2000 epoch for the equatorial
coordinate system. We express the natural logarithm as $\ln(x)$ and the
base-10 logarithm as $\log(x)$.

##  2 Low-redshift sample of remnants

We compile a sample of remnant radio galaxies from a range of large-sky radio
surveys to form the observational basis for the later analysis in this work
(section 2.1). The radio galaxies in this parent sample are cross-matched with
infrared identified host galaxies, which are then matched to obtain optical
passband photometries and redshifts, and to infer host galaxy properties
(section 2.2). We then use literature-established methods to select remnant
lobes from our parent radio galaxy sample (section 2.3). Finally, following
the method of Quici et al. (2022), we measure the observable attributes of
each remnant lobe needed to constrain its energetics (section 2.4).

###  2.1 Radio galaxy selection

We outline the steps followed to select a parent sample of radio galaxies from
existing large-sky surveys. We impose the selection criteria that the radio
sources are brighter than $0.5$ Jy at 150 MHz (i.e.,
$S_{150\,\mathrm{MHz}}\geqslant 0.5$ Jy) and larger than four arcminutes in
size ($\theta>4^{\prime}$); this ensures the measured radio source attributes
are not heavily degraded either by (relatively) poor surface brightness
sensitivity or a lack of spatial sampling across the lobe.

####  2.1.1 Radio survey availability

We consider sources selected by both the GaLactic and Extragalactic All-sky
MWA (GLEAM; Hurley-Walker et al., 2017) survey and the Rapid ASKAP Continuum
Survey (RACS; McConnell et al., 2020) to ensure radio images are present at
two or more frequencies, whilst maximising the potential sample of remnant
radio galaxies. These two surveys both cover the following declination range:
$-85^{\circ}\leqslant\mathrm{Decl.}\leqslant 30^{\circ}$. We supplement these
observations with images from other large-sky surveys with partial coverage of
this region: the TIFR GMRT Sky Survey (TGSS; Intema et al., 2017), the NRAO
VLA Sky Survey (NVSS; Condon et al., 1998) and the Very Large Array Sky Survey
(VLASS; Lacy et al., 2019). The properties of these surveys are summarised in
Table 1. Importantly, the observational parameter space provided by the GLEAM
survey, principally the low-frequency coverage and low spatial-resolution,
minimises the deselection of aged remnant lobes.

Survey | Frequency | Sensitivity | Resolution | Decl. limit
---|---|---|---|---
| (MHz) | (mJy/beam) | (′′) |
GLEAM | 72-216 | $\sim 7$ | $120$ | $\leqslant+30^{\circ}$
TGSS | 151 | $\sim 5$ | $25$ | $\geqslant-53^{\circ}$
RACS-low | 887 | $\sim 0.25$ | $\sim 25$ | $\leqslant+30^{\circ}$
NVSS | 1400 | $\sim 0.45$ | $45$ | $\geqslant-40^{\circ}$
VLASS | 3000 | $\sim 0.12$ | $3$ | $\geqslant-40^{\circ}$
Table 1: Properties of the large-sky radio surveys used in this work. Each row
represents a unique radio survey (column 1), for which the observing
frequencies (column 2), root mean square sensitivities (column 3), angular
resolutions (column 4), and declination limits (column 5) are shown.

####  2.1.2 Integrated flux density cut: $S_{150\,\mathrm{MHz}}\geqslant 0.5$
Jy

We begin the parent sample compilation by considering all radio components
brighter than $S_{150\,\mathrm{MHz}}\geqslant 0.5$ Jy. We conduct this step
using the GLEAM survey due primarily to its enhanced surface brightness
sensitivity to large-scale structure. As a preliminary step, we incorporate
the GLEAM 4 Jansky catalogue (G4Jy; White et al., 2020a, b), which offers a
complete GLEAM-selected sample of radio galaxies brighter than 4 Jy at 151
MHz, and also offers robust host galaxy associations. Importantly, all sources
in the G4Jy catalogue meet the integrated flux criteria of our sample, and are
thus filtered through to our candidate sample. To avoid any duplicates of
these sources, we remove components comprising the radio galaxies present in
the G4Jy sample from the GLEAM catalogue prior to further source selection.
This is done using a cone search centered at the coordinates of the G4Jy
centroid, setting the radius of the search equivalent to reported angular
size.

We now need to identify all GLEAM-selected radio sources with an integrated
flux density of $0.5\leqslant S_{150\,\mathrm{MHz}}(\mathrm{Jy})\leqslant 4$.
For each component in the GLEAM catalogue, we use the fitted $72-216$ MHz
spectral index ($\alpha_{72}^{216}$) and reference 200 MHz integrated flux
density ($S_{200\,\mathrm{fit}}$) to compute the 150 MHz integrated flux
density; we use this value over that measured from the $154$ MHz GLEAM wide-
band, in order to reduce the measurement uncertainty due to white noise.
Components with flux densities of $S_{150\,\mathrm{MHz}}\geqslant 0.5$ Jy are
added to our candidate sample. However, recognising that extended
($\theta>4^{\prime}$) radio galaxies can appear as multiple components even in
the GLEAM catalogue (e.g., see Fig. 1), we must also consider radio sources
for which the sum of their component flux densities exceeds our required
threshold. To do this, we perform a cone search within the remaining GLEAM
catalogue to internally cross match all components within a
$\theta=18^{\prime}$ radius; i.e., based on the largest angular size reported
in the G4Jy catalogue. Component groups with a total flux density
$S_{150\,\mathrm{MHz}}\geqslant 0.5$ Jy are also added to our candidate
sample.

Following this approach, we identify a total of 3071 candidate radio sources
brighter with $S_{150}\geqslant~{}0.5$ Jy.

####  2.1.3 Extended radio galaxy cut: $\theta\geqslant 4^{\prime}$

(a) J1528+0544

(b) J1317-2422

Figure 1: Example overlays used to manually inspect the GLEAM-selected radio
sources. Grey-scale background images represent the AllWISE W1 bands. Each
radio source is shown at 200 MHz as seen by GLEAM (solid black contours), and
as seen by RACS-low at 887 MHz (semi-transparent images); the shape of their
restoring beams are shown in the lower-left corner by the black and purple
ellipses, respectively. For Fig. 1(a), a close-up view around the host galaxy
coordinates, as seen at 3 GHz by VLASS, is shown in the upper-right corner.
For Fig. 1(b), the coordinates of the host galaxy are marked by the cyan
cross.

We next need to identify which of our radio source candidates are: (1)
extended radio galaxies larger than $\theta\geqslant 4^{\prime}$, and (2) are
not comprised of components associated with unrelated radio sources.
Considering the size of this sample is just on the limit of what can be
handled manually, we perform this step largely through a visual inspection to
ensure our classifications are robust. However, we perform two initial cuts to
remove compact radio galaxies from our sample, as we now briefly discuss.

We use the Tool for OPerations on Catalogues And Tables (TOPCAT) software,
with a search radius of $\theta=2^{\prime}$ (i.e., the size of the GLEAM
resolution at 200 MHz), to perform an internal cross match within the TGSS
catalogue (for $\mathrm{Decl.}\geqslant-40^{\circ}$) to search for all
isolated components within this radius. We then filter the TGSS components for
sources with a integrated-to-peak flux density ratio less than $1.3$; i.e.,
definitely compact. These TGSS objects are cross-matched with the GLEAM
catalogue, and we remove the associated GLEAM components from the sample. We
perform the same internal cross match within the NVSS catalogue to search for
isolated components in this survey, and cross-match these outputs with the
GLEAM catalogue. By calculating a two-point spectral index between $151\,$MHz
and $1.4$ GHz, GLEAM components are excluded if the measured spectral index is
greater than zero (i.e., brighter at higher frequencies). This is a
conservative cut, considering that radio galaxy lobes should have a spectral
index at least steeper than $\alpha=-0.5$, and is designed to exclude compact
objects such as gigahertz peaked spectrum objects and blazars.

The remaining radio galaxy candidates are visually inspected to check that the
components are associated with an actual radio galaxy and meet the angular
size cut ($\theta\geqslant 4^{\prime}$). We use the coordinates of their radio
centroid to create $0.5^{\circ}\times 0.5^{\circ}$ postage stamps from the
GLEAM, TGSS, RACS-low and VLASS radio images. For multi-component source
candidates, we take the radio centroid as the mean of the individual component
coordinates. Importantly, the (relative) high resolution and surface
brightness sensitivity of the RACS-low radio imaging make it possible to
inspect the radio morphologies of each candidate radio source in sufficient
detail. We use the SAOImageDS9 (Joye & Mandel, 2003) software to manually
measure the largest angular size (LAS) across the lobes. Sources are discarded
if the measured LAS fails to reach the required angular size criterion.
Objects are also discarded if their GLEAM component(s) are associated with
unrelated sources. This is true for both single-component sources, i.e., where
the singular GLEAM component is in fact multiple radio sources, as well as
multi-component sources, i.e., where the GLEAM components themselves are
revealed to be unrelated. These decisions are aided with infrared images
indicating potential host galaxies locations (discussed in section 2.2), and
also aided with the discarding of diffuse radio emission associated with
nearby star-forming galaxies (SFGs).

Finally, following this approach, we arrive at a parent sample of 795 radio
galaxies with $S_{151}>0.5$ Jy and $\theta\geqslant 4^{\prime}$. To prepare
our sample for further analysis, we compile the integrated flux densities for
each source from the five GLEAM wide-bands, RACS-low and NVSS (for
$\mathrm{Decl.}\geqslant-40^{\circ}$).

###  2.2 Local sub-halo environment

The properties of the environment within which the source expands need to be
quantified to successfully model the evolutionary history of a radio galaxy.
Specifically, we require estimates for the host galaxy redshift and stellar
mass, as we briefly discuss. The RAiSE code uses semi-analytic galaxy
evolution model outputs to identify plausible dark matter halo properties
based on the measured stellar mass and cosmological evolution at the observed
redshift. The median relationship between these variables is used to provide
point estimates for the halo mass in this work (see Turner, 2017). Meanwhile,
the shape of the ambient gas density profile in RAiSE is informed by X-ray
observations scaled to the properties of the chosen dark matter halo (Turner &
Shabala, 2015; Turner et al., 2023).

In this section, we cross match our radio galaxy sample to the Wide-field
Infrared Survey Explorer (WISE) all-sky survey (section 2.2.1), and use the
WISE object IDs to automatically cross match our sample against various online
data bases in order to obtain redshift measurements (section 2.2.2) and
optical passband photometries to calculate value-added products (section
2.2.3). Unless otherwise stated, all database cross matching outlined in this
section is performed using the VizieR catalogue
service111https://vizier.cds.unistra.fr/.

####  2.2.1 The AllWISE Data Release

We perform our initial host galaxy cross match with respect to the AllWISE
Data Release (Cutri et al., 2021) publicly available on
VizieR222https://vizier.cds.unistra.fr/viz-bin/VizieR?-source=II/328. We
choose this survey considering it covers the full declination range occupied
by our sample, and also probes a wavelength range well-suited to the detection
of dust-obscured host galaxies. Additionally, its most sensitive passband (W1,
$\lambda=3.4\,\mu$m) has a high angular resolution of $\sim 4^{\prime\prime}$,
comparable to that of VLASS. We create similar postage stamps to those
discussed in section 2.1.3 for the $\lambda=3.4\,\mu$m band of the WISE survey
(Wright et al., 2010). These images were used to identify the coordinates of
the likely host galaxy based on the orientation of the radio emission with
respect to the galactic environment. For radio sources with a clear radio core
and/or central jet, we took the host galaxy as that which coincided with the
origin of the radio emission. We used the highest resolution radio data
available, which at minimum was provided by RACS-low, however in many cases
was improved on by VLASS imaging. For radio sources without a clear radio
core, we recorded the likely coordinate of the host galaxy based on the
orientation of the diffuse lobe plasma with respect to the surrounding
galaxies. We note, although the host galaxy identification is often
challenging for remnant radio galaxies (see e.g., Brienza et al., 2017; Quici
et al., 2021), we did not encounter any cases where the host galaxy associated
was ambiguous. We obtain the WISE object ID for the candidate host galaxies
identified using the WISE W1 images (see Fig. 1).

####  2.2.2 Host galaxy redshifts

Redshift information is obtained for our radio sources using the astroquery
(Ginsburg et al., 2019) python module to cross match our host galaxy
coordinates with the NASA/IPAC Extragalactic
Database333https://ned.ipac.caltech.edu/ (NED). For each radio source, we
perform a $\theta=2^{\prime}$ cone search around the coordinates of the host
galaxy, and filter the outputs for objects with catalogued redshifts in NED.
Through this approach, redshifts are found for 540 out of the original 795
host galaxies. We also record the classification of each redshift type,
whether the redshifts are spectroscopic or photometric, based on the original
references for each catalogued redshift. Overall, our sample was broken down
into 329 spectroscopic, and 211 photometric redshifts. Our redshift
distribution is shown in Fig. 2. As a sample cut, we limit the redshift
distribution to $0.02<z<0.2$, to reduce the contamination from sources with
photometric redshifts (but see also discussion in section 3.1.2). At
$z\geqslant 0.2$ our photometric redshifts become more numerous than
spectroscopic which, considering the redshift is needed to perform the
parameter inversion, may introduce an unwanted systematic into the fitted
energetics.

Figure 2: The redshift, stellar and halo mass distributions for our parent
radio galaxy sample. The two scatter plots show the derived halo masses (upper
panel) and stellar masses (lower panel) as a function of redshift.
Spectroscopic and photometric redshift are denoted using blue and orange
colors, respectively. Each histogram shows the one-dimensional distribution in
the parameter adjacent to the plot. The grey bars denote the full sample
(spectroscopic and photometric).

Since these redshifts are heterogeneously-selected, we can not comment on the
redshift completeness of our sample. However, we clarify that the method
undertaken for this work is not reliant upon volume-complete radio galaxy
samples. In section 4, the predictions made by forward modelling a given jet
power and lifetime model are made specific to the intrinsic property
distributions of the selected remnant sample. This includes the redshift
distribution, which in principle means that the results should not be biased
by redshift incompleteness.

####  2.2.3 Host galaxy stellar masses

We derive the stellar mass of each host galaxy based on the empirical
relationship between the rest-frame $(g-i)$ colour and mass-to-light ratio
($M_{*}/L_{i}$) (Taylor et al., 2011). The stellar mass is given by their
linear best fit (e.g., see their equation 8):

$\log[M_{*}/M_{\odot}]=-0.68+0.70(g-i)-0.4m_{i}-2\log\Big{(}\frac{d_{\mathrm{L}}}{10\,\mathrm{pc}}\Big{)},$
(1)

where $m_{i}$ is the apparent magnitude in the photometric $i$-band (AB
system), $d_{\mathrm{L}}$ is the luminosity distance, and the $1\sigma$
uncertainty corresponding to the derived $M_{*}/L_{i}$ is $\sim 0.1$ dex.

However, no single optical survey provides coverage of our entire radio galaxy
sample. As such, we resort to compiling measurements from various surveys,
specifically the Panoramic Survey Telescope and Rapid Response System (Pan-
STARRS1; Chambers et al., 2016), the Sloan Digital Sky Survey Data Release 12
(SDSS DR12; Alam et al., 2015), and the Skymapper Southern Sky Survey (SMSS;
Wolf et al., 2018). Using the coordinates of the WISE object IDs obtained
above, we perform a cone search of each database using a
$\theta=5^{\prime\prime}$ search radius to retrieve the best positional match
for each object. For each survey, we obtain the AB magnitudes in the $g$ and
$i$ passbands. For regions below a declination of $-40^{\circ}$, photometries
are taken from SMSS. Photometries are otherwise taken from PanSTARRS, however
for a small fraction of these objects ($\sim 1\%$), the PanSTARRS photometry
is either missing in one of the passbands, or has a poorer signal to noise
ratio than that in SDSS. In these two cases, optical magnitudes are sourced
from SDSS instead.

In order to minimise any potential systematic biases that could lead to
inconsistent stellar-mass estimates from multiple surveys for the same galaxy,
we apply corrections to the SMSS and SDSS $g$ and $i$ magnitudes to align them
with PanSTARRS, which serves as the reference due to its overlap with both
surveys. The galaxies in our sample that appear in the overlapping region with
PanSTARRS are used for independent linear corrections of the $g$ and $i$ bands
in each survey.

Figure 3: Consistency check for stellar masses derived using the $(g-i)$
colors from the PanSTARRS, SMSS and SDSS surveys. In each panel the PanSTARRS-
based mass estimate (horizontal axis) is compared to the SMSS-based mass
estimate (left panel) and the SDSS-based mass estimate (right panel). The
standard deviation of the differences calculated between the two mass
estimates is shown in the label for each panel, and represented by the grey-
shaded region.

In Fig. 3, we compare the stellar mass estimates (equation 1) of each survey
for sources in the overlapping region with PanSTARRS. For the comparison of
both the SMSS and SDSS surveys against PanSTARRS, we find the linear best fit
to the data and compute the standard deviation of the residual. Reassuringly,
our results demonstrate good agreement between all three stellar mass
estimates. While the scatter is clearly larger for the SMSS-PanSTARRS
comparison, both comparisons show a near 1:1 relationship, and thus do not
appear to be systematically offset. We assume the same holds true for the
entire sample. To account for the noisier SMSS data, we assign a $0.32$ dex
uncertainty (e.g., the $1\sigma$ deviation of the scatter) on the SMSS-derived
stellar masses.

Our stellar mass distribution, and corresponding halo mass distribution
derived from semi-analytic galaxy evolution model outputs (Turner & Shabala,
2015), is shown in Fig. 2. The vast majority of galaxies have stellar masses
between $10^{10.5}$ and $10^{12}$ M⊙, however, a few have extreme values
fitted. The stellar mass of these galaxies is likely poorly constrained with
two wavelength bands, however, as the mass is solely used as a prior to inform
the shape of the host galaxy ambient medium, we simply assume semi-plausible
values (either $10^{10.5}$ or $10^{12}$ M⊙) for these objects.

###  2.3 Classifying remnant radio galaxies

We classify remnant radio galaxies based on the findings of Brienza et al.
(2017), Mahatma et al. (2018), Jurlin et al. (2020) and Quici et al. (2021),
who show that remnants can display a wide range of observed properties. In
section 4, we will apply forward modelling to directly constrain the
energetics of the sample of remnant radio galaxies compiled here. For this
reason, provided our methods do not preferentially select specific sub-
populations, we prioritise our remnant selection for reliability (i.e., how
confident we are in our classification) rather than completeness (i.e., how
certain we are we have identified all the remnants in our parent sample).
Furthermore, despite the observational differences between remnant and
restarted radio galaxies, the latter are equally valuable in this work.
Restarted radio galaxies are defined as those in which multiple episodes of
activity are observed concurrently. Importantly, the associated remnant lobes
are no different to those associated with an inactive nucleus, and thus can be
modelled through a RAiSE-based parameter inversion provided the two epochs of
radio activity do not overlap spatially. With this in mind, we use the cutouts
created in section 2.1.3 and adopt several criteria to perform our remnant
(and restarted) radio galaxy selection.

Prior to our remnant selection, for each source we calculate several “flags”
of remnant classification based on its observed radio properties. First, radio
sources are identified as remnants if they do not exhibit radio emission from
the core. We classify a core detection if radio emission at the coordinates of
the host galaxy is significant above the $3\sigma$ level, where $\sigma$
denotes the local root mean square (rms) noise. This step is performed
uniformly across our sample using the RACS-low data, but additionally includes
VLASS data for $\text{Decl.}\geqslant-40^{\circ}$ sources. We classify the
core status as ambiguous if radio emission is present above this level, but is
diffuse and/or amorphous in morphology (discussed in the following paragraph).
Additionally, using the VLASS and/or RACS-low data as above, sources are
identified as remnants if they do not show compact features such as jets or
hotspots. Next, we calculate a ‘low surface brightness’ classification metric
based on the mean 151 MHz lobe surface brightness (in mJy/arcmin-2). Using the
RACS-low data, we use the
polygon_flux444https://github.com/nhurleywalker/polygon-flux code (Hurley-
Walker et al., 2019) to trace the $3\sigma$ footprint of each radio source,
and thus quantify the area projected by the lobes. Together with the 151 MHz
GLEAM integrated flux density, we then evaluate an average surface brightness,
and identify sources that fall below $50\,$mJy/arcmin-2 as remnants. We
calculate an ‘ultra-steep spectrum’ classification metric, based on a global
spectral index fitted to the lobe spectrum. For each source, we fit the global
spectral energy distribution (between the lowest and highest frequency
available) in log-log space using a simple linear model, and identify remnants
as those sources with a fitted spectral index steeper than $\alpha=-1.2$.

We are now able to select remnant radio galaxies from our parent sample using
the classification scheme described above. The ‘absent radio core criterion’
(Quici et al., 2021), in principle, can classify sources as remnant if they do
not exhibit significant radio emission at the host galaxy. However, Jurlin et
al. (2021) showed that faint radio emission from the core, perhaps due to the
jets not completely shutting off or star formation at lower redshifts, can
still coincide with remnant radio galaxy lobes. For sources where the core is
classified as ‘ambiguous’, we adapt the ‘morphological’ criterion (Brienza et
al., 2017) and classify sources as remnant if the following two conditions are
met: (1) their lobes are ‘flagged’ as having ‘low surface brightness’, and (2)
the sources do not exhibit any jets/hotpots. In support of this step, we
select the known remnant radio galaxy associated with NGC 1534 (Hurley-Walker
et al., 2015; Duchesne & Johnston-Hollitt, 2019), which has unambiguously
remnant lobes, yet also shows weak radio emission at the host galaxy. Finally,
in line with the restarted radio galaxy classification techniques of Jurlin et
al. (2020), we additionally select remnant lobes if a clear radio core is
observed, however their lobes exhibit the following properties: their lobes
are classified as ‘low surface brightness’ and as ‘ultra-steep spectrum’, and
do not exhibit any compact features (jets/hotspots).

Following these methods, we obtain a final sample of 79 remnant radio galaxies
within $0.02<z<0.2$.

###  2.4 Measuring the radio source attributes

We need to measure a set of attributes from the observed radio images of our
selected remnant radio galaxies in order to perform a Bayesian parameter
inversion for their energetics (section 5). In the RAiSE framework, these
observed radio source attributes are compared to those extracted from
synthetic radio sources using a maximum likelihood algorithm. It follows,
then, that the quality and number of measured radio source attributes (unique
constraints) must be adequate to constrain the underlying parameters. For this
work, we implement the result of Quici et al. (2022), who showed that
attributes related to the observed surface brightness distribution of remnant
lobes help to constrain the duration of the remnant phase.

With this in mind, and given the quality of the radio data available to us
(see section 2.1.1), the following radio source attributes are of particular
interest to us:

  * •

a reference radio luminosity;

  * •

a largest linear size;

  * •

the two extent attributes, which describe the lobe surface brightness
distribution (e.g., Quici et al., 2022);

  * •

a low frequency spectral index; and

  * •

a measure of spectral curvature.

We define the reference radio luminosity, $L_{\nu}$, as the 154 MHz integrated
flux density measured for each remnant. The decision to measure this attribute
using GLEAM data is motivated by the sensitivity to low surface brightness
emission and relatively low calibration uncertainty (e.g., approximately 8%;
Hurley-Walker et al., 2017). However, for three of our selected remnants,
blending with unrelated components resulted in an uncertain measurement of
their $S_{154}$. For these sources, we measured their reference radio
luminosity at 887 MHz using RACS-low data (which can be equally well predicted
by RAiSE).

To measure the largest linear size and the two extent attributes, we use the
method of Quici et al. (2022), who parameterise these attributes by fitting
the lobe surface brightness profiles with skewed Gaussian functions (e.g., see
their section 4.2.2). As part of their analysis, they demonstrated that the
extent attributes are mostly affected by the spectral ageing of the oldest
plasma, meaning that higher-frequency measurements are more constraining. We
therefore measure the extent attributes at 887 MHz, considering that the RACS-
low observations provide the optimal trade-off between frequency, sensitivity,
and spatial resolution; for example, while NVSS probes higher frequencies, the
restoring beam is also larger, and thus does not offer much constraining power
here.

The low-frequency spectral index, $\alpha_{\mathrm{low}}$, is measured using
the five GLEAM 30 MHz-wide bands (ranging between 88 MHz and 200 MHz). In the
ideal case, the spectral index is measured across all five bands. However, for
a handful of cases, blending from unrelated components (which becomes
increasingly worse towards lower frequencies) makes this measurement
problematic. For two such sources, the low-frequency spectral index was
measured between 154 MHz and 200 MHz. To derive the spectral index, spectral
data points were fit in log-log space (flux density and frequency) using a
simple linear fit.

Finally, to obtain the spectral curvature, SPC, we measured a higher-frequency
two-point spectral index, $\alpha_{\mathrm{high}}$, and calculated the
difference with $\alpha_{\mathrm{low}}$ (see Quici et al., 2022). To measure
this attribute, we used one of two ways. If the radio source was detected by
NVSS, $\alpha_{\mathrm{high}}$ was calculated between 887 MHz and 1.4 GHz.
Consequently, if the radio source fell below the NVSS declination limit
($\text{Decl. = }-40^{\circ}$), or above the declination limit but was
undetected by NVSS, we calculated $\alpha_{\mathrm{high}}$ between $200$ MHz
and $887$ MHz. We found that all sources detected by NVSS were also detected
by RACS-low, which is unsurprising considering the deeper and comparatively
lower-frequency observations provided by RACS-low.

With this approach, the radio source attributes are measured for each remnant
in our sample. Using these attributes, we now seek to model the evolutionary
histories of our remnants, and thus constrain their energetics. This is
discussed in the following section.

##  3 Constraining the sample energetics

The energetics of an observed radio source can be constrained via a RAiSE-
based parameter inversion – a method that exploits the ability of the RAiSE
model to convert from intrinsic parameter space into observed attribute space.
In the following section, we describe the adoption of this method to quantify
the energetics for the sample of remnants compiled in the previous section.
Briefly summarised, synthetic radio sources are generated across a multi-
dimensional grid of input parameters that are expected to influence the
observed properties (section 3.1). The parameter space is inverted by
comparing synthetic attributes from the simulated surface brightness maps to
their observed counterparts, whilst carefully modelling the redshift and
survey limitations of each radio source. (section 3.2). We use these results
to explore the observed distributions in the fitted parameters (section 3.3).

###  3.1 Creating the RAiSE simulation grid

The RAiSE simulations that underlie this work play a key role in both
estimating the energetics of each observed remnant (section 3.2.2), and
generating mock remnant populations required later in this work (section 4).
We assume these remnants have a Fanaroff & Riley 1974 Type-II morphology in
their active phase, given this class is associated with higher jet powers and
thus more likely to be detected after the cessation of the jet (Turner, 2018).
Towards both of these aims, we must ensure that the input parameter space
covers the full breadth of plausible values, and that the grid resolution in
each dimension is at least comparable to the associated fitting uncertainty.
Below, a detailed setup of the RAiSE simulation is outlined.

####  3.1.1 Initial setup

To predict the evolutionary histories of a radio source, we use the Turner et
al. (2023) ‘jet+lobe’ RAiSE model, for which the associated python code is
publicly available on GitHub555github.com/rossjturner/RAiSEHD. Radio sources
are simulated using the RAiSE_run function, which generates spatially-resolved
images of the lobe synchrotron radiation for an set of model inputs. Two such
inputs, frequency and resolution, control the properties of the output images.
For computational efficiency, we use resolution=poor, which sets the number of
Lagrangian particles used in the surface brightness calculation equal to
1,792,000. To match the combined frequency coverage provided by the GLEAM,
RACS-low, and NVSS radio surveys, the synchrotron radiation is synthesised at
88, 151, 200, 887 and 1400 MHz. The inputs related to the dynamics and
energetics of the lobes represent dynamically relevant parameters that we will
need to fit for each source and are discussed in the following section. The
initial axis ratio could arguably be varied, but due to computational
constraints, we assume a fixed value of $A_{0}=2.5$ consistent with
observations of 3CRR Fanaroff & Riley (1974) Type-II radio galaxies (Turner et
al., 2018b); crucially, Turner et al. (2018a, their section 6.5) show that
this parameter has a very weak influence on the fitted energetics. Other input
parameters for the jets/lobes, e.g., those describing the internal structure
of the jet, are set to default values that have been calibrated against
hydrodynamical simulations (see Table 1 of Turner et al., 2023). Finally, in
these simulations we underlying the radio source environments based on a halo
mass, which places the radio source at the centre of a spherically-symmetric,
double-beta gas density profile (Turner & Shabala, 2015).

####  3.1.2 The input parameter space

The RAiSE-based parameter inversion uses attributes measured for an observed
radio source to infer its intrinsic parameters. In this approach, a hyper-
dimensional grid of order $k$ is required to fit $k$ parameters, where each
dimension is a list of uniformly-gridded values (in either log- or linear-
space). Our main goal is to fit the energetics of each radio source, meaning
we need to fit the jet power and active age. However, we also want to account
for other intrinsic properties that will have an influence on the observed
radio source properties, thus ensuring that the fitted energetics are not
systematically incorrect. These properties are: the remnant ratio (fractional
time spent in a remnant phase;
$R_{\mathrm{rem}}=t_{\mathrm{rem}}/(t_{\mathrm{on}}+t_{\mathrm{rem}})$), halo
mass, equipartition factor of the lobes (measure of the magnetic field
strength), and the energy injection index. The intrinsic parameter space is
thus created by uniformly gridding each aforementioned parameter between an
upper/lower limit; see Table 2 for details. This results in a $k=6$
dimensional grid of mock radio sources.

Table 2: The properties of the parameter space upon which the hyper-dimensional RAiSE simulation grid is created. Each row represents one of the six dimensions of the simulation grid, which is gridded between an upper/lower limit (column 3) using a fixed resolution (column 5), and comprises a total of $n$ values (column 6). This parameter space results in a total of $8.63\times 10^{7}$ synthetic radio sources.  Input parameter | Symbol | Value(s) | Unit | Resolution | $n$
---|---|---|---|---|---
Jet power | $Q$ |  $10^{36.5}$ \- $10^{40}$ | W |  $0.1$ dex | $36$
Active age | $t_{\mathrm{on}}$ |  $10^{6.5}$ \- $10^{8.5}$ | Myr |  $0.05$ dex | $41$
Remnant ratio | $R_{\mathrm{rem}}$ |  $0$ \- $0.7$ | - | $0.025$ | $29$
Halo mass | $M_{\mathrm{halo}}$ |  $10^{11.8}$ \- $10^{15}$ | $\text{M}_{\odot}$ |  $0.2$ dex | $16$
Equipartition factor | $q$ |  $10^{-2.7}$ \- $10^{-0.7}$ | - |  $0.1$ dex | $21$
Injection index | $s$ |  $2.2$ \- $2.7$ | - | 0.1 | $6$

Ideally, to fit the intrinsic parameters for our remnant sample, a unique
RAiSE simulation grid would be created at the redshift of each remnant.
However, the RAiSE model outputs are only a weak function of redshift:
notably, radiative losses from the inverse-Compton upscattering of CMB photons
depends on redshift as $(1~{}+~{}z)^{4}$. Considering that the redshifts of
our observed remnant sample is low ($z<0.2$) we can reasonably approximate
this loss mechanism with a simulation at a single redshift, $z=0.1$. This
reduces the dimensionality of our parameter space thereby reducing the total
computation time by a factor of 79 (our sample size). The observable
attributes are thus calculated for the observed redshift of each remnant by
scaling the RAiSE simulation outputs (run at $z=0.1$) prior to convolution
with the relevant survey beam.

Following this approach, the RAiSE simulation grid is constructed by
simulating the radio source evolutionary histories for each point in the input
($k=6$ dimensional) parameter space. Each point contains spatially-resolved
maps of the synchrotron radiation expected at the frequencies specified in
section 3.1.1. The integrated radio luminosities (for each input frequency)
and the true physical size (measured from hotspot to hotspot) are also are
calculated by default for each RAiSE model.

###  3.2 Modelling the intrinsic parameters

The energetics of each remnant can now be estimated via a parameter inversion
using the RAiSE model as the theoretical expectation. To do this, the
synthetic images must be scaled to match the redshift of the remnant, then be
degraded to match the sensitivity and angular resolution of the radio
observations; the radio source attributes are then measured for each synthetic
image (section 3.2.1). The observed and synthetic radio source attributes are
compared for each synthetic image using a maximum likelihood algorithm and
weighted using Bayesian priors for properties of the host environment and lobe
magnetic field. This gives probability density functions for the $k=6$
intrinsic parameters of each remnant in our sample, e.g., jet power, active
age and off-time (section 3.2.2).

####  3.2.1 Synthetic radio source attributes

The RAiSE simulations are run at a single redshift for computational
efficiency, as discussed in section 3.1.2. These two-dimensional (2D)
synchrotron-emission maps are recast into the frame of the observed radio
source based on the measured redshift; in this process, we convert the pixel
values of the output image from a luminosity (W/Hz) into a flux density (Jy),
and similarly convert the scale of the pixel from a physical size (kpc) into
an apparent size (′′). The scaled outputs are convolved with a 2D Gaussian
kernel with a scale and orientation matching the angular resolution of the
relevant survey (e.g., to mimic a RACS-low observation, images at 887 MHz are
convolved with a circular beam of diameter $25^{\prime\prime}$). These
convolved outputs are converted into synthetic surface brightness images
(pixel values in units of Jy/beam) by dividing through by the area of the
restoring beam. We note that the synthetic images do not include a random
noise component as our sample is selected to have a good signal to noise
ratio. Simulated radio source attributes, matching those measured for each
remnant (section 2.4), are extracted from the synthetic, survey-degraded
surface brightness images (see also Quici et al., 2022). Following this
approach, we obtain a set of attributes for the mock radio source generated by
each combination of potential intrinsic radio galaxy properties – these sets
of attributes can be compared to those of the observed remnants in a parameter
inversion.

####  3.2.2 Bayesian parameter estimation

The energetics of each remnant are fitted using a maximum likelihood algorithm
based on a comparison of the observed radio source attributes (section 2.4) to
the ($k=6$) hyper-dimensional grid of attributes for each combination of
intrinsic parameters modelled using RAiSE (section 3.2.1). Following Quici et
al. (2022), the Akaike information criterion (AIC) corresponding to each RAiSE
model is calculated by comparing the observed and synthetic attributes via the
likelihood function (e.g., their equation 9). The optimal model (or set of
intrinsic parameters) is selected as that which minimises the AIC; the
probability of each model is calculated relative to this best-fit. These
probabilities only encode the agreement between simulation and observation,
and do not reflect the Bayesian priors that exist for several parameters. For
this work, the evaluated probabilities are multiplied through by two priors:
the logarithmic equipartition factor, $\log q$, which is known to be
approximately Gaussian-distributed ($\mu\approx-1.73$, $\sigma\approx 0.53$;
Turner et al., 2018b); and the halo mass, for which we obtained a measure of
the mean and uncertainty in section 2.2 based on the $(g-i)$ stellar mass.

Following this approach, the energetics of each remnant in our sample are
constrained within the hyper-dimensional parameter space. The expectation and
uncertainty in each of the $k=6$ parameters are derived from the marginal
probability density functions (PDFs); i.e., the probabilities obtained sum
over all but the parameter of interest. We take the 50th percentile of the PDF
as our estimate of the best fit, and the 16 and 84th percentiles as our
estimates of the 1$\sigma$ confidence interval. We examine the fitted
properties of our remnant sample in the following section.

###  3.3 Observed jet-power and lifetime distributions

We have now constrained the energetics of each observed remnant using the
RAiSE-based parameter inversion, notably the jet kinetic power and duration of
the active phase. We use Monte Carlo (MC) statistics to characterise the shape
and uncertainty in the fitted jet power and lifetime distributions for our
sample. Each MC realisation consists of $N=79$ jet power and active age
estimates based on randomly sampling the probability density functions for
each remnant. The resulting jet power and lifetime distributions for each
realisation are then binned and normalised (area equal to unity) in a two-
dimensional parameter space (i.e., in jet power and active age), where the
bins are based on the gridded parameter space of the RAiSE simulation grid
(Table 2). We perform 1000 MC realisations to obtain an estimate for the
expectation and uncertainty in the simulated counts for each bin in the jet
power–active age, $(Q,t_{\mathrm{on}})$, parameter space. Our results are
shown in Fig. 4.

Figure 4: The observed marginal PDFs in jet power (lower right) and active age
(upper left) fitted for our remnant sample. Monte Carlo simulations are used
to build up the observed sample statistics in each parameter (grey shading).
The two-dimensional grid (lower right) shows the observed joint PDF in
$(Q,t_{\mathrm{on}})$-space. The colour-map indicates the mean probability
density, $p$, of sources at each pixel.

Importantly, due to the sample selection bias, these distributions do not
represent the true underlying statistics for these parameters. We can see that
the observed jet power distribution peaks at $Q\approx 10^{38.5}$ W, and shows
an observed dearth of low-powered sources; this turn-over is likely due to
Malmquist bias (see also section 4 below) – the preferential selection of
bright objects in sensitivity-limited surveys. A similar effect, primarily due
to the minimum angular size cut, can be seen for the active age distribution,
which peaks at $t_{\mathrm{on}}\approx 10^{8.25}$ yrs.

In section 4, we will attempt to decode this sample selection bias by
comparing mock remnant population with known jet power and lifetime functions
to our observed distribution. To this end, it is worth noting that there is no
obvious correlation between the jet power and active age in Fig. 4, suggesting
we can reasonably ignore higher-order interactions between these parameters in
our work; i.e., we can assume the jet power and active age are only weakly
correlated. Regardless, in the following section, we will investigate their
distributions in $(Q,t_{\mathrm{on}})$-space, rather than the marginal PDFs,
to ensure we constrain these parameters independently of each other.

##  4 Simulating the observed energetics

The jet power and lifetime distributions fitted for our remnant sample in
section 3.3 will not represent the underlying population statistics due to the
sample selection biases. We outline a process to resolve this problem by
comparing the energetics of a mock population generated with known underlying
distributions – with sample selection biases applied – to those fitted for our
observed remnant sample. First, we generate mock remnant catalogues for an
assumed jet power and lifetime distribution, and model how the selection bias
modifies the observed energetics (section 4.1). Second, we must control for
the confounding variables (i.e., redshift, halo mass, equipartition and
injection index) so that the predictions made for a given jet power and
lifetime distribution can be meaningfully compared to our actual sample
(section 4.2).

###  4.1 Modelling the sample selection bias

Mock catalogues of remnant radio galaxies are created by weighting the
probability of the parameter space investigated in the RAiSE simulation. To
model the effect selection bias has on the observed energetics, two separate
weighting schemes are required. First, the simulation needs to be weighted by
an underlying jet power and active age distribution (section 4.1.1); we refer
to this weighting as $W_{\mathrm{jet\text{-}age}}$, which controls the assumed
shape of the underlying distributions we will later constrain. Second, to
compute the selection bias, $W_{\mathrm{selection}}$, a binary weighting is
assigned to each synthetic source depending on whether or not the source would
pass the sample selection criteria (section 4.1.2).

####  4.1.1 The input jet power and lifetime distributions

We begin with the assumption that a simple power-law model is sufficient to
describe the distributions of the radio-loud AGN jet powers and lifetimes.
Observational constraints on the AGN radio luminosity function (RLF) show that
the shape of the RLF is a broken power law (e.g., see Sabater et al., 2019;
Franzen et al., 2021). Considering that, to first order, the radio luminosity
is correlated with jet power (noting that the age and environment also have a
confounding impact, e.g., Godfrey & Shabala 2013; Hardcastle & Krause 2013),
and that our sample preferentially selects high-powered remnants (see Fig. 4),
we assume our assumption holds for the jet power distribution. Meanwhile,
simulations of feedback-regulated black hole accretion show that the
Eddington-scaled mass accretion rate,
$\dot{m}=\dot{M}_{\mathrm{BH}}/M_{\mathrm{Edd}}$, follows a power spectrum
consistent with pink noise (e.g., see Novak et al., 2011; Gaspari et al.,
2017). These simulations predict that the jets triggered at the black hole are
described by a lifetime function consistent with a $p(t_{\mathrm{on}})\propto
t_{\mathrm{on}}^{-1}$ power-law slope666This function describes the likelihood
of a jet being triggered with an active age $t_{\mathrm{on}}$, not the
relative probability of observing such a source.. Based on our discussion in
section 3.3, we assume that the two distributions (jet power and active age)
are seeded independently; that is, that no intrinsic correlation exists
between these two quantities. Future work based on a larger sample could
investigate more complex model assumptions including correlations with the
host galaxy mass (see section 6).

Based on these arguments, we model the underlying jet power and lifetime
distributions as the following probability density functions (PDFs):

$\displaystyle p(Q)\,dQ=\left({Q}/{Q_{0}}\right)^{a}\,dQ,$ (2a) $\displaystyle
p(t_{\mathrm{on}})\,dt_{\mathrm{on}}=\left({t_{\mathrm{on}}}/{t_{0}}\right)^{b}\,dt_{\mathrm{on}},$
(2b)

where $Q_{0}$ and $t_{0}$ are arbitrary constants that ensure the total
probability is unity between the lower and upper bounds (see Table 2). The
coefficients $a$ and $b$ give the power-law indices for the jet power and
lifetime distributions of each jet outburst (i.e., not considering the
likelihood of observing).

However, as the input parameter space for the RAiSE simulation is uniformly
gridded in $\log Q$ and $\log t_{\mathrm{on}}$, weighting the simulation using
the above expressions would fail to account for the non-uniform grid spacing
when transformed to linear space. This is important considering that, for
example, the grid bin at $\log(t_{\mathrm{on}}/\text{yrs})=7.5$ covers is an
order of magnitude greater range of ages than that at
$\log(t_{\mathrm{on}}/\text{yrs})=6.5$. As such, the above expressions are
adapted by the bin correction as follows:

$\displaystyle p(\log Q)\,d\log Q=\left(Q/Q_{0}\right)^{a+1}Q_{0}\ln 10\,d\log
Q,$ (3a) $\displaystyle p(\log t_{\mathrm{on}})\,d\log
t_{\mathrm{on}}=\left(t_{\mathrm{on}}/t_{0}\right)^{b+1}t_{0}\ln 10\,d\log
t_{\mathrm{on}}.$ (3b)

We will not provide results in the form of these equations – these are only
necessary to weight the simulation outputs to the correct, temporally uniform
sampling assumed in equations 2a and 2b.

The remnant radio lobes modelled in this work are observed exclusively after
the jet has switched off (by definition)777Restarted radio sources included in
our sample may have an active jet confined to the host galaxy – in this case,
our methods are explicitly modelling the outburst associated with the remnant
lobe.. The active age is correctly weighted for uniform sampling in the
linearly-spaced active age parameter (equation 2b), but we have not considered
if this holds for the remnant phase, parameterised through the remnant ratio
$R_{\mathrm{rem}}=t_{\mathrm{rem}}/(t_{\mathrm{on}}+t_{\mathrm{rem}})$. We
find the following relationship between a uniform increment in the linearly-
spaced remnant ratio (see Table 2) and a uniform increment in off-time:

$dt_{\mathrm{rem}}=\frac{(t_{\mathrm{on}}+t_{\mathrm{rem}})^{2}}{t_{\mathrm{on}}}\;dR_{\mathrm{rem}}.$
(4)

This equation implies that simulation outputs for a given remnant ratio must
be weighted by the factor
$(t_{\mathrm{on}}+t_{\mathrm{rem}})^{2}/t_{\mathrm{on}}$ to ensure uniform
sampling in the off-time. We apply this correction to our simulation outputs.

The weightings in equations 3a, 3b and 4 specify the weighting required both
to encode the underlying jet power and lifetime functions, and to convert from
the sampling in our simulation grid to a temporal uniform sampling. The
function $W_{\mathrm{jet\text{-}age}}$ represents the weights calculated in
this manner for each point in the $k=6$ parameter space, enabling us to
generate mock populations, albeit without selection effects applied as yet.

####  4.1.2 Applying the selection criteria

The selection biases inherent to the observed sample (section 2.1) are
replicated for each mock source in the RAiSE simulation grid by assessing
their synthetic surface brightness images against our sample selection
criteria. Three key criteria merit consideration: (1) the minimum integrated
flux density; (2) the minimum surface brightness; and (3) the minimum angular
size measured along the radio lobes. The modelling of each selection criteria
is discussed below, and is intended to replicate the steps outlined in section
2 for the observed sample.

To construct the parent catalogue of radio sources, the GLEAM catalogue was
used to select radio sources brighter than $0.5$ Jy at 150 MHz (i.e.,
$S_{150}\geqslant 0.5\,$Jy). For each synthetic radio source, the 150 MHz
integrated flux density is calculated based on the forward-modelled radio
luminosity and the given mock redshift. As such, we assume that the synthetic
sources fainter than this threshold would not have made it into our sample.
Considering that the GLEAM catalogue contains only the components with a
signal to noise larger than $5\sigma$, where $\sigma$ is the GLEAM local rms
in Jy/beam, this sets the condition for the minimum surface brightness
criteria. To do this, the synthetic 150 MHz synchrotron emission maps are
convolved to match the GLEAM resolution at 150 MHz. Using these surface
brightness maps, we assume that synthetic sources with a peak flux density
fainter than the $5\sigma$ GLEAM sensitivity would not have been detected, and
thus will be missed by our selection. Finally, our sample contains radio
sources with a projected angular size greater than $\theta=4^{\prime}$ in
RACS-low. To model this selection criteria, the synthetic 887 MHz synchrotron
emission maps are reprojected for a given mock redshift, and convolved to
match the RACS-low resolution at 887 MHz – we assume all the observed radio
sources lie in the plane of the sky. Considering only the regions of the lobes
brighter than $750\,\mu$Jy/beam ($\approx 3\sigma$), we use the surface
brightness maps to discard sources whose largest angular size falls below our
sample threshold. To model the bias due to these sample selection criteria, we
construct a binary selection grid as follows:

$W_{\mathrm{selection}}=\bigg{\\{}\\!\begin{array}[]{ll}1,&\text{selected}\\\
0,&\text{not selected}\end{array},$ (5)

where the synthetic source is weighted either as unity (the source passes the
selection criteria), or zero (the source fails one or more of the selection
criteria).

We note, the framework above only models the selection due to the parent
sample selection criteria (i.e., section 2.1), and does not consider the
potential selection biases due to the remnant selection methods we have used.
For this work, we assume that as soon as a source switches off it will show
observable signatures consistent with a remnant radio galaxy; in other words,
the remnant radio galaxies forward modelled by our mock catalogues are ones
that we should be able to classify. This assumption may not hold for extremely
young remnants, for example those in which ${t_{\mathrm{rem}}}$ is smaller
than the light travel time along the jet,
${t_{\mathrm{light\text{-}travel}}}$888For a distance of $100$ kpc,
information travelling at a speed $0.3c$ would take $\sim 1$ Myr to reach the
hotspot from the SMBH., however, we do not expect this to be a problem for the
large off-times probed by our sample.

Following this approach, the selection of each source is evaluated for the
entire RAiSE simulation. Importantly, by evaluating
$W_{\mathrm{jet\text{-}age}}$ and $W_{\mathrm{selection}}$, predictions for
the observed jet power and lifetime distributions can be made for any
underlying assumptions and sample selection criteria. To make predictions for
the observed energetics, the RAiSE simulation grid is sampled using the
element-wise multiplication of the two weighting schemes: i.e.,
$W_{\mathrm{jet\text{-}age}}\circ W_{\mathrm{selection}}$.

This approach is demonstrated in Fig. 5 assuming a mock population of $z=0.1$
remnant galaxies for jet power and active age distributions of $p(Q)\propto
Q^{-1.2}$ and $p(t_{\mathrm{on}})\propto t_{\mathrm{on}}^{-0.5}$ respectively.
The sample selection bias against lower-powered and shorter-lived sources is
clearly seen: these sources are fainter and have a smaller physical sizes, and
are thus precisely the objects that would be preferentially deselected from
our sample. Higher-powered jets produce stronger lobe magnetic fields, in turn
resulting in shorter fading timescales during the remnant phase (e.g., Turner
et al., 2018b; Turner, 2018; Quici et al., 2022); this effect is likely
responsible for the observed flattening towards $Q=10^{40}$ W. Similarly, the
downturn observed towards the longest-lived sources is likely a combination of
two factors: the expansion of the source driving down the lobe surface
brightness (e.g., Turner, 2018; Shabala et al., 2020), and the decrease in the
integrated radio luminosity due to Rayleigh-Taylor mixing of the lobe plasma
with the environment (e.g., see Fig. 3 of Turner & Shabala 2015). Importantly,
these results show us sample selection biases can greatly transform the
observed nature of the true jet power and active age distributions.

Our approach so far, however, has neglected the impact the remaining sample
properties will have on these results. This is discussed in the following
section.

###  4.2 Controlling for the global properties of the sample

The range of plausibly-selected jet power and lifetime distributions are
strongly confounded by the redshift, in addition to the halo mass,
equipartition factor, and injection index. The selection biases computed for
our mock catalogues must therefore take into account how these parameters are
distributed for our sample. To do this, we initially apply a zeroth-order
correction: assume that their underlying distributions are the same as what we
observe (section 4.2.1). However, after applying the sample selection biases
the PDFs for each parameter will likely differ from their input distributions;
we therefore apply higher-order corrections by iteratively perturbing the
underlying distributions of these parameters until – upon applying the
selection effects – they match the observed sample (section 4.2.2).

####  4.2.1 Weighting the initial priors

We use the observed statistics in each parameter to weight their corresponding
mock distributions; e.g., the observed redshift distribution acts as a prior
to weight the redshift distribution of the synthetic sources, as we now
describe in detail. The observed redshifts are binned to to the closest
matching mock redshift, along a uniformly spaced axis from $(0,0.2]$ using a
resolution of $dz=0.005$999The RAiSE simulations run at $z=0.1$ are scaled to
these mock redshifts for computational efficiency, as discussed in section
3.1.2 for our remnant sample fitting.. This assignment of the observed
redshifts is repeated using 1000 Monte Carlo (MC) realisations of the each
redshift within its measurement uncertainties. The resulting distribution is
normalised to unity; this normalised PDF for the redshift is taken as the
prior to weight our simulation outputs. We note, empty bins are assigned a
weighting of zero, meaning there is a null probability of finding a source
with that particular parameter value.

The observed distributions of the other confounding variables are built up
using MC statistics as for the redshift, with the simulation grid for each
parameter specified in Table 2. The uncertainties on these confounding
variables are taken as the fitted error when modelling the evolutionary
histories of each remnant in section 3.2.2; i.e., unlike the redshift, we do
not have robust measurements for these parameters before the RAiSE-based
inversion. In this way, we have characterised the observed distributions of
each confounding variable and converted these to prior PDFs for each axis of
the RAiSE simulation grid.

Figure 5: The underlying distributions (grey bars) in the jet power (left
panel) and active age (right panel) with $p(Q)\propto Q^{-1.2}$ and
$p(t_{\mathrm{on}})\propto t_{\mathrm{on}}^{-0.5}$ respectively. This mock
catalogue is simulated at $z=0.1$, $\log(H/\text{M}_{\odot})=14$, $\log
q=-1.7$ and $s=2.2$. We filter this catalogue through our sample selection
criteria (orange; section 2.1). The prior weightings for the redshift, halo
mass, and other confounding variables, are derived from the sample
distributions as either a zeroth-order correction (light blue; section 4.2.1)
or are updated iteratively to match the sample distributions, acting as a
first-order correction (dark blue; section 4.2.2).

The zeroth-order correction to our mock population is simply to weight each
axis of the RAiSE simulation grid by the observed distributions in these
confounding variables. We apply these prior PDFs to the example population
considered in the previous section (shown in Fig. 5), demonstrating the effect
of considering the remaining properties of our sample. The most notable
difference is the turn-off towards the highest-powered jets; a feature that is
also observed for our remnant sample (Fig. 4). This is likely due to some
higher-powered jets not forming lobes in the lowest-density environments
(i.e., balistic expansion; Falle, 1991; Turner & Shabala, 2023).

####  4.2.2 Controlling for confounding variables

We have so far applied a zeroth-order correction to our confounding variables
by assuming their underlying distributions match those observed for our
sample. However, this is unlikely to be true as, like for the jet power and
active age, the distributions in these parameters will be moderated by the
sample selection effects. We demonstrate this problem as follows: we generate
two mock catalogues with different underlying energetics, weight their
confounding variable PDFs following the method outlined in the previous
section, and finally apply the selection criteria of our sample. The mock
redshift and halo mass PDFs – resulting after applying the selection biases –
are compared to the observed distributions in Fig. 6, showing the selection
bias does indeed moderate the shape of the observed distributions. The output
redshift and halo mass PDFs do not match their input distributions, implying
that the observed statistics of our remnant sample must have been selected
from a currently-unknown underlying distribution.

Figure 6: The observed redshift (left panel) and halo mass (right panel)
distributions predicted for two separate mock catalogues (solid red and dot-
dashed orange steps) created for different jet power and active age
distributions (see exponents in figure label). Mock catalogues are filtered by
the selection criteria of our sample, and are weighted by the distributions
observed for the redshift and the confounding variables (grey bars). The
redshift bins plotted here are larger than those of the mock simulation grid
to aid this visualisation.

The manner by which the outputs differ to the input is consistent with our
expectation. The output redshift PDF is skewed towards lower-redshifts: for a
fixed jet power and lifetime distribution, the highest redshift bin will have
the narrowest range of plausibly-selected evolutionary histories. As for the
halo mass distribution, the bias is evident in lower mass environments where
high-powered jets are less likely to form lobes. In Fig. 6, it is also
apparent that the selection bias acting on the redshift and halo mass is
dependent upon the underlying energetics (red versus orange lines). The model
with more low-powered and short-lived sources ($a=-1.4$ and $b=-1.0$) tends to
be deselected at higher redshifts (fewer bright sources) and in higher mass
haloes (fewer jets meeting the angular size cut). Overall, these findings
imply that the underlying distributions for the confounding variables must
have higher-order corrections applied such that the resulting survey-selected
PDFs match those of our observed sample.

The higher-order corrections are derived by iteratively applying corrections
to the confounding variable PDFs as in section 4.2.1. We apply the zeroth-
order correction as before, then, for each confounding variable, we compare
the error in each bin (Table 2) between the observed and mock PDFs for a
proposed jet power and lifetime distribution; i.e.,
$N_{\mathrm{observed}}/N_{\mathrm{mock}}$. This ratio quantifies how
significant the selection bias is at each bin, which varies somewhat smoothly
across the parameter space. The prior PDFs are updated by the inverse of this
ratio to up-weight values that were under-represented after the previous
correction. The mock catalogues are re-sampled with the updated prior PDFs for
each of the confounded variables. This process can be repeated over several
iterations, however, we find that a single round of (first-order) corrections
is sufficient to bring the mock distributions in agreement with our remnant
sample.

Finally, with these corrections, we are able to produce mock catalogues that
are correctly weighted for their redshift and other confounding-variable
distributions. The predictions made by these catalogues, e.g., see Fig. 5, can
now compared to the observed jet power and active age distribution of our
sample, in order to constrain their underlying distributions. This is explored
in the following section.

##  5 Results and discussion

We now can separate the impact of selection bias from our observed jet power
and lifetime distributions using the method outlined in section 4. This allows
us to infer their underlying distributions (section 5.1) and compare these
findings with previous work in the literature (section 5.2). In section 5.3,
we discuss our results in the context of jet triggering and fuelling
mechanisms, and then we suggest techniques to apply our approach to the active
radio source population in section 5.4.

###  5.1 Constraining the underlying distributions

Following the methodology described in section 4, we generate mock catalogs in
a two-dimensional grid of proposed jet power and lifetime distributions, where
each cell corresponds to the power-law indices $(a,b)$. The power-law index
for the jet power distribution is binned into discrete intervals in
$-2.6\leqslant a\leqslant-0.6$ using a precision of $\delta a=0.1$. Similarly,
the power-law index for the lifetime distribution is partitioned into
intervals in $-2\leqslant b\leqslant-0.5$ with a precision of $\delta b=0.05$.

We first investigate whether our mock catalogues (for the proposed of
underlying distributions) exhibit discernible differences when influenced by
the selection biases. We compare the predicted jet power distribution between
two scenarios: a steep jet power distribution with $a=-2.0$ (many low-powered
jets) and a shallow distribution with $a=-1.0$ (fewer low-powered jets), and
consider all age models within $-2\leqslant b\leqslant-1$. Similarly, we
compare the predicted lifetime distribution between the following two
scenarios: a steep distribution with $b=-2.0$ (many short-lived jets) and a
shallower distribution with $b=-1.0$ (fewer short-lived jets), and consider
all jet power models within $-2\leqslant a\leqslant-1$. Fig. 7 presents our
findings, which are compared with the observed sample statistics from Fig. 4.
Encouragingly, despite the stringent selection criteria, the shape of the
predicted distributions (after selection biases are applied) are measurably
different: those corresponding to shallower underlying distributions exhibit a
turn-over at approximately one order of magnitude greater than the steeper
underlying distribution.

Figure 7: Predictions for the observed distributions in jet power (left panel)
and active age (right panel) made by various mock populations. The left panel
shows predictions made by a steep ($a=-2.0$, red) and shallower ($a=-1.0$,
blue) jet power model, evaluated for all age distributions. The right panel
shows predictions made by a steep ($b=-1.0$, brown) and shallower ($b=0$,
teal) ageing model. evaluated for all jet power models with
$1\leqslant|a|\leqslant 2$. The mock distributions are smoothed using a $0.3$
dex kernel.

Fig. 7 further illustrates the intertwined relationship between the underlying
jet power and active age distributions, and their selection-biased
counterparts. For a given jet power model, the shape of the observed jet power
distribution is somewhat influenced by the underlying jet lifetime model, and
vice versa. This indicates that while the observed distribution primarily
constrains its corresponding underlying distribution (i.e., observed jet
powers constrain the jet power model and observed active ages constrain the
lifetime model), there should in principle be weaker constraints stemming from
the complementary distribution (i.e., observed active ages aiding in
constraining the jet power model, and vice versa).

To constrain the evolutionary jet model, the mock jet power and active age
distributions are binned in $(Q,t_{\mathrm{on}})$-space to match the binning
of the observed distributions (see Fig. 4). The likelihood function is
evaluated by considering the mock and observed constraints at each point in
this parameter space, from which we calculate the AIC for two free parameters.
This calculation assumes that each data point (pixel) offers an independent
measurement on the number of sources in that particular $(Q,t_{\mathrm{on}})$
bin, which is not strictly true considering that there is a correlation
between the number of sources in neighbouring bins. We note however that the
uncertainty on the observed $(Q,t_{\mathrm{on}})$ distribution is estimated
via Monte Carlo simulation (see section 3.3), and encodes the effect of this
correlation. We select the optimal jet model as that which minimises the AIC,
and we compute the relative probability of each model as
$p_{i}=e^{(\mathrm{AIC}_{i}-\mathrm{AIC}_{\mathrm{best}})/2}$. The 2D
probability density function (PDF) is mapped across $(a,b)$-space (Fig. 8);
the PDF shows some covariance, implying that the jet power and lifetime
distribution are not fitted completely independent of each other.

Figure 8: The relative fitting probabilities evaluated for each jet power and
lifetime model. The corresponding probability density function (PDF) mapped in
$(a,b)$-space is shown in the bottom left panel. A sequential colour-map is
used to indicate the probability for each model, and is shown using a
logarithmic stretch. Contours denote the $1\sigma$, $2\sigma$ and $3\sigma$
confidence intervals. The adjacent panels demonstrate the PDF marginalised for
$a$ (upper left panel), and $b$ (lower right panel), as shown by the solid
black curves. The red dot-dashed curves show the Gaussian approximations of
the marginal PDFs.

The parameter estimates are quantified using their 1D marginal PDFs, shown in
Fig. 8. We approximate the marginal PDFs using a Gaussian function to
characterise the most probable value and the $1\sigma$ standard deviation;
this reduces the need to sample the parameter space more finely at large
computational expense. As such, we constrain the power-law slopes of the
underlying jet power and lifetime functions as follows:

$\displaystyle p(Q)\propto Q^{-1.49\pm 0.07},$ (6a) $\displaystyle
p(t_{\mathrm{on}})\propto t_{\mathrm{on}}^{-0.97\pm 0.12}.$ (6b)

The resulting survey-selected jet power and lifetime distributions are
compared to the observed distributions in Fig. 9. We can see that the shape of
the mock active age distribution agrees with observation; the mock
distribution follows the same initial rise in number counts as the observed
distribution towards longer-lived sources, and exhibits a similar turnover at
$t_{\mathrm{on}}=10^{8.25}$ yrs. Within uncertainties, the general shape of
the mock jet power distribution also appears consistent with observation; the
number density of sources at the higher and lower-powered ends are in
agreement with each other, and the turn overs in each distributions are
comparable. However, several potential inconsistencies should be pointed out.
Firstly, the turnover in the mock jet power distribution ($Q=10^{39}$ W)
occurs roughly an order of magnitude higher than that in the observed
distribution ($Q\approx 10^{38.2}$ W). Secondly, the rate at which the mock
distribution falls off towards low-powered sources is not as steep as what
appears to be the case for our observed sample. In fact, by manually
inspecting the predictions made our mock catalogues (but see Fig. 7), we find
that no model is able to simultaneously replicate each of these features.
Flatter values of $a$ give the correct fall-off towards low $Q$, but turn over
at much higher jet powers. Steeper values of $a$ predict a closer turn over to
what is observed, but over predict the number density of low-powered objects.

Figure 9: Same as Fig. 4, but overlaying the predictions made by our best-fit
jet power and lifetime model; i.e., $p(Q)\propto Q^{-1.5}$ and
$p(t_{\mathrm{on}})\propto t_{\mathrm{on}}^{-1.0}$.

We speculate that one potential explanation could be that the true jet power
function flattens below a certain jet power, $Q_{*}$, which would be expected
if the broken power-law nature of the AGN radio luminosity function is driven
directly by shape of the jet power function. Such a model should in principle
provide a solution to the two inconsistencies outlined above; by flatting
below some $Q_{*}$, the expected fall-off towards low-powered sources should
be sharper than what is currently shown by our mock catalogue in Fig. 4.
Considering our relatively low sample size $(N=79)$, we have chosen not to fit
a broken power-law jet power model for this work. Nevertheless, under the
assumption that the jet power function is indeed more appropriately described
by a broken power law, our results here show a promising future application
for radio galaxy populations revealed by next-generation radio surveys.

###  5.2 Comparison with the literature

####  5.2.1 AGN radio luminosity functions

Historically, our strongest constraints towards the shape of the underlying
jet power function have come from the radio-loud AGN radio luminosity function
(RLF). Extensive work exists within the literature that constrains the AGN
RLF, and we focus on several recent results to examine the consistency with
our model expectations.

Using the GLEAM radio survey, Franzen et al. (2021) measured the 200 MHz RLF
for local (median $z\approx 0.064$) radio galaxies. By fitting their observed
AGN RLF using a broken power-law model, they constrain a power-law slope of
$\alpha_{\mathrm{RLF}}=-1.76$ above a critical radio luminosity
($L_{\star}=10^{25.76}$ W Hz-1). Using a theory-driven first-order
approximation for the conversion between radio luminosity and jet power,
$Q\propto L_{\nu}^{6/7}$ (Willott et al., 1999; Kaiser & Best, 2007), this
suggests that the power-law slope in their corresponding jet power function is
$a\approx-2$. Though considerably steeper than the power-law index constrained
by our method ($a=-1.49\pm 0.07$), comparing these values must be exercised
with caution. Since the observed luminosity of a radio source is strongly
confounded by its jet power, age, and ambient environment (Hardcastle &
Krause, 2013; Shabala & Godfrey, 2013; Yates et al., 2018), the resulting
scatter in the AGN RLF prevents the direct inference of the jet power
function. The effects of a selection bias may also contribute towards the
apparent mismatch between the fitted jet power slopes. Although powerful radio
sources will have bright radio luminosities, such sources will also have large
physical sizes, which may contribute to a deselection from brightness-limited
radio surveys.

Intriguingly, the 150 MHz AGN RLF measured by Shimwell et al. (2019) – who
used radio sources selected as part of the first data release of the LOFAR
Two-Metre Sky Survey (LoTSS DR1; Shimwell et al., 2019) – appears flatter at
the high-luminosity end as compared to Franzen et al. (2021). It is entirely
possible that this is driven by their low-number statistics at brighter
luminosity bins, since their sample probes a smaller cosmological volume at
low redshifts. However, taking their result at face value would imply a
flatter jet power distribution at high jet powers, and thus qualitatively in
better agreement with our fitted value.

####  5.2.2 Literature forward modelling studies

Forward modelling has recently been employed in two published works to
constrain the bulk properties of radio-loud AGNs based on the distributions in
their observed attributes, as we discuss below.

Hardcastle et al. (2019) used the dynamical model of Hardcastle (2018) to
apply forward modelling to radio-loud AGNs unveiled by LoTSS DR1. As part of
their analysis, they tested the validity of two different lifetime models in
which the observed active ages were distributed uniformly in either linear age
($b=0$ in this work), or in logarithmic age ($b=-1$ in this work). Their jet
power function was assumed to be $p(Q)\propto Q^{-1}$. To examine the validity
of their lifetime models, the linear size distributions were forward modelled
for each model assumption, and compared to the observed size distributions for
three radio luminosity bins. Their main finding (see their section 4.2) is
that neither lifetime function was able to simultaneously explain the observed
size distributions across all luminosity bins. Their linear age model gave
excellent agreement with the size distributions in the brightest two bins,
however vastly under-predicted the fraction of small ($<100\,$kpc) radio
sources in the faintest bin. On the other hand, at the highest two luminosity
bins, their log-uniform lifetime function vastly over-predicted the compact
fraction, and simultaneously under-predicted the observed size distributions.

Shabala et al. (2020) employed RAiSE forward-modelling to constrain the jet
power and lifetime functions for a LOFAR-selected sample of radio galaxies.
For a given jet power and lifetime function, they forward-modelled radio
galaxy populations, applied sample selection biases, and compared their
predicted angular size and flux density distributions to those observed for
their sample. In particular, they tested power-law and constant-age models,
and found that both could plausibly explain the observed properties of their
active radio sources. To break down the degeneracy between these two lifetime
functions, they compared the mock remnant and restarted fractions with
observation, and found that power-law lifetime function offered much better
agreement. In particular, they found that a model with $p(Q)\propto Q^{-1}$
and $p(t_{\mathrm{on}})\propto t_{\mathrm{on}}^{-1}$ offered excellent
agreement with their observed sample, with $p(t_{\mathrm{on}})\propto
t_{\mathrm{on}}^{-0.5}$ also offering good agreement.

We now examine the consistency between our results and those of the two
aforementioned studies, but note that differences are expected due to
differences in the methodologies we have adopted and the radio galaxy samples
we have used. We constrain the same optimal lifetime model as Shabala et al.
(2020), i.e., $p(t_{\mathrm{on}})\propto t_{\mathrm{on}}^{-1}$, though our
results show their $p(t_{\mathrm{on}})\propto t_{\mathrm{on}}^{-0.5}$ model is
incompatible with the observed active age distribution. These findings are
also tentatively in agreement with Hardcastle et al. (2019). Considering that
we have directly measured the full lifetime of the jets in these sources, we
would expect the fitted active ages to occupy a single value (within fitting
uncertainties) if their underlying lifetimes were that of a constant age. On
the other hand, both Hardcastle et al. (2019) and Shabala et al. (2020) find
considerably flatter jet power models, which may be a result of the sample
selection. In section 5, we speculated that the true jet power function may be
more appropriately described by a broken power law, such that a flattening
would occur below some critical value. The jet powers fitted by Shabala et al.
(2020) for their sample show a peak in their distribution at $Q\approx
10^{37.7}$ W (see their fig. 3), almost an entire order of magnitude lower
than our sample. We therefore speculate that the flatter jet power
distribution constrained in their work could simply be a consequence of their
sample probing a lower-powered portion of the jet power function.

###  5.3 Jet triggering and fuelling mechanisms

A key result of this work is the confirmation that the jet lifetime function
is consistent with $p(t_{\mathrm{on}})\propto t_{\mathrm{on}}^{-1}$, implying
that most radio sources are short-lived. The active age distribution fitted
for our remnant sample offers direct evidence against constant age models
($>$$5\sigma$ confidence level), and is in support of the conclusions of
Shabala et al. (2020). This finding begs the question: what jet triggering
mechanisms lead to more short-duration outbursts? Simulations of feedback-
regulated black hole accretion consistently predict a time-varying mass
accretion rate. The associated power spectrum is consistent with pink noise,
corresponding to a lifetime function $p(t_{\mathrm{on}})\propto
t_{\mathrm{on}}^{-1}$ (e.g., see Novak et al., 2011; Gaspari et al., 2017).
This agreement with our jet lifetime function, and earlier work by Shabala et
al. (2020), strongly suggests the jet active age distribution is moderated by
feedback-regulated accretion.

We note that numerous authors have suggested that different jet triggering
mechanisms may mediate the lifetime functions for different radio source
populations (e.g., Pimbblet et al., 2013; Kaviraj et al., 2015; Marshall et
al., 2018). However, our remnant sample probes a population with jet kinetic
powers approximately an order of magnitude higher than that of Shabala et al.
(2020); the pink noise spectrum indicative of feedback-regulated accretion is
encoded over several orders of magnitude in AGN jet power. Regardless, neither
sample considers significant populations of lower-powered Fanaroff & Riley
(1974) Type-I sources ($Q\lesssim 10^{37}$ W) so our conclusions may be
limited to the higher-powered Type-II class.

###  5.4 Application to active radio galaxy populations

While the approach taken for this work focused solely on remnant radio
galaxies, our methods for simulating mock radio source populations can be
generalised to the wider radio-loud AGN population; specifically to active
radio galaxies, which make up $\sim 80\%$ of the observed radio source
population (e.g., see Jurlin et al., 2020). As mentioned in section 1, the age
measured for an active radio galaxy represents the age at which the source is
observed, $\tau$, not the true duration of its active phase. This pushed
Turner & Shabala (2015) to assume that, on average, active radio sources will
be observed halfway through their active lifetime, i.e.,
$\tau=t_{\mathrm{on}}/2$. However, an an alternative approach to this problem
could involve directly simulating an observed $\tau$ distribution. For a mock
radio source, $\tau$ can be simulated by uniformly sampling, in linear space,
the range bounded by $t=0$ and $t=t_{\mathrm{on}}$. For any underlying jet
lifetime function, a corresponding $\tau$ distribution can thus be simulated
and filtered by the selection criteria, as outlined in section 4. The
simulated, selection-biased age distribution can be directly compared to the
source ages, $\tau$, fitted for active radio galaxies. Of course, a major
caveat with this method is that a large observational sample (e.g.,
$\sim$$1000$ sources) would be required to confidently average over the random
noise resulting from sampling each active source at some unknown point in its
active phase – the statistical ensemble resulting from numerous objects is
needed to inform the active age in this case.

##  6 Conclusion

We have presented an approach to constrain the jet power and active age
functions for active galactic nuclei (AGN) jet outbursts by decoding
observations of remnant radio galaxies. These objects uniquely measure the
full duration of the active phase, thus placing firm constraints on the jet
lifetime function. We compiled a sample of remnant radio galaxies using large-
sky multi-frequency radio data with a fixed selection criteria, i.e.,
$S_{150}\geqslant 0.5$ Jy and $\theta\geqslant 4^{\prime}$, and identified
remnant lobes using variety of literature-established methods (section 2). The
energetics of each source were constrained using the RAiSE dynamical model;
simulated radio source attributes were compared with observation through a
parameter inversion, which allowed fitting of, amongst other parameters, their
jet power and active age (section 3). However, due to sample selection biases,
the fitted jet power and active age distributions do not represent their true
underlying population statistics. To combat this, we simulated mock remnant
populations for assumed (power law) models of their jet powers and lifetimes,
and filtered by the sample selection criteria (section 4). Mock catalogues of
remnant populations were generated for a range of jet power and lifetime
functions and compared to the observed jet powers and active age
distributions. We constrained the optimal pairing of jet kinetic power and
lifetime functions and compared our results to the theoretical expectations
for jet triggering and fuelling mechanisms.(section 5). We summarise our
findings as follows:

  * •

We constrain the jet kinetic power and lifetime functions as
$p(Q)\propto~{}Q^{-1.49\pm 0.07}$ and
$p(t_{\mathrm{on}})\propto~{}t_{\mathrm{on}}^{-0.97\pm 0.12}$, respectively.

  * •

Our results are consistent with feedback-regulated accretion simulations,
which predict a power-law spectrum in black hole accretion rates equivalent to
a lifetime function of $p(t_{\mathrm{on}})\propto t_{\mathrm{on}}^{-1}$.

  * •

Our jet power function is broadly consistent with the high-luminosity end of
the AGN radio luminosity function (Shimwell et al., 2019; Franzen et al.,
2021). We speculate that the true jet power function may be better described
by a broken power-law model.

  * •

Our jet lifetime function aligns well with the findings of Hardcastle et al.
(2019) and Shabala et al. (2020). However, our approach allows us to impose
more precise constraints on the range of plausible lifetime functions, notably
ruling out their $p(t_{\mathrm{on}})\propto~{}t_{\mathrm{on}}^{-0.5}$ model.

Despite the focus of this work being on remnant radio galaxies, our
methodology is applicable to active radio galaxy populations (section 5.4).
Considering these objects dominate flux and volume limited samples, their
statistical significance is an attractive property that should be exploited in
similar analyses. In the advent of new generation radio surveys such as LoTSS,
the Evolutionary Map of the Universe (EMU; Norris, 2011), and eventually those
produced by the Square Kilometer Array, constraints on the global energetics
of radio-loud AGNs will dramatically improve. Not only will these surveys
provide large, statistically-significant samples of radio galaxies, their
sensitivity and angular resolution will allow for much lessened selection
biases. This will allow, for example, the testing of higher-order models
(e.g., double power-law jet kinetic power models), and investigating whether
the correlations exist between the jet kinetic powers and lifetimes, or with
properties of the host galaxy. Furthermore, with the inclusion of higher-
quality optical/infrared (e.g., the LOFAR-WEAVE survey; Smith et al., 2016)
and X-ray (e.g., eROSITA; Predehl et al., 2021) data, our methodology can be
extended to explore the environmental impact on AGN fuelling.

## Acknowledgements

BQ acknowledges a Doctoral Scholarship and an Australian Gov- ernment Research
Training Programme scholarship administered through Curtin University of
Western Australia. NHW is supported by an Australian Research Council Future
Fellowship (project number FT190100231) funded by the Australian Government.
This work was supported by resources provided by the Pawsey Supercomputing
Research Centre (https://ror.org/04f2f0537) with funding from the Australian
Government and the Government of Western Australia.

## Data Availability

The authors confirm that the data supporting the findings of this study are
available in cited external data archives or presented in the article.
Processed data products underlying this article will be shared on reasonable
request to the authors. The relevant code is publicly available as disclosed
when referenced in the paper.

## References

  * Alam et al. (2015) Alam S., et al., 2015, ApJS, 219, 12
  * Alexander & Hickox (2012) Alexander D. M., Hickox R. C., 2012, New Astronomy Review, 56, 93
  * Best et al. (2005) Best P. N., Kauffmann G., Heckman T. M., Brinchmann J., Charlot S., Ivezić Ž., White S. D. M., 2005, MNRAS, 362, 25
  * Bower et al. (2006) Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645
  * Brienza et al. (2016) Brienza M., et al., 2016, A&A, 585, A29
  * Brienza et al. (2017) Brienza M., et al., 2017, A&A, 606, A98
  * Chambers et al. (2016) Chambers K. C., et al., 2016, arXiv e-prints, p. arXiv:1612.05560
  * Churazov et al. (2002) Churazov E., Sunyaev R., Forman W., Böhringer H., 2002, MNRAS, 332, 729
  * Condon et al. (1998) Condon J. J., Cotton W. D., Greisen E. W., Yin Q. F., Perley R. A., Taylor G. B., Broderick J. J., 1998, AJ, 115, 1693
  * Crockett et al. (2012) Crockett R. M., et al., 2012, MNRAS, 421, 1603
  * Croft et al. (2006) Croft S., et al., 2006, ApJ, 647, 1040
  * Croton et al. (2006) Croton D. J., et al., 2006, MNRAS, 365, 11
  * Cutri et al. (2021) Cutri R. M., et al., 2021, VizieR Online Data Catalog, p. II/328
  * Daly (2009) Daly R. A., 2009, ApJ, 696, L32
  * Daly (2016) Daly R. A., 2016, MNRAS, 458, L24
  * Daly (2021) Daly R. A., 2021, MNRAS, 500, 215
  * Dasyra & Combes (2011) Dasyra K. M., Combes F., 2011, A&A, 533, L10
  * Dasyra & Combes (2012) Dasyra K. M., Combes F., 2012, A&A, 541, L7
  * Dubois et al. (2014) Dubois Y., Volonteri M., Silk J., 2014, MNRAS, 440, 1590
  * Duchesne & Johnston-Hollitt (2019) Duchesne S. W., Johnston-Hollitt M., 2019, Publ. Astron. Soc. Australia, 36, e016
  * Duchesne et al. (2020) Duchesne S. W., Johnston-Hollitt M., Zhu Z., Wayth R. B., Line J. L. B., 2020, Publ. Astron. Soc. Australia, 37, e037
  * Dugan et al. (2017) Dugan Z., Gaibler V., Silk J., 2017, ApJ, 844, 37
  * Emonts et al. (2014) Emonts B. H. C., Piqueras-López J., Colina L., Arribas S., Villar-Martín M., Pereira-Santaella M., Garcia-Burillo S., Alonso-Herrero A., 2014, A&A, 572, A40
  * English et al. (2016) English W., Hardcastle M. J., Krause M. G. H., 2016, MNRAS, 461, 2025
  * English et al. (2019) English W., Hardcastle M. J., Krause M. G. H., 2019, MNRAS, 490, 5807
  * Fabian (2012) Fabian A. C., 2012, ARA&A, 50, 455
  * Fabian et al. (2003) Fabian A. C., Sanders J. S., Allen S. W., Crawford C. S., Iwasawa K., Johnstone R. M., Schmidt R. W., Taylor G. B., 2003, MNRAS, 344, L43
  * Falle (1991) Falle S. A. E. G., 1991, MNRAS, 250, 581
  * Fanaroff & Riley (1974) Fanaroff B. L., Riley J. M., 1974, MNRAS, 167, 31P
  * Forman et al. (2005) Forman W., et al., 2005, ApJ, 635, 894
  * Franzen et al. (2021) Franzen T. M. O., et al., 2021, Publ. Astron. Soc. Australia, 38, e041
  * Gabor & Bournaud (2013) Gabor J. M., Bournaud F., 2013, MNRAS, 434, 606
  * Gaspari et al. (2017) Gaspari M., Temi P., Brighenti F., 2017, MNRAS, 466, 677
  * Ginsburg et al. (2019) Ginsburg A., et al., 2019, AJ, 157, 98
  * Godfrey & Shabala (2013) Godfrey L. E. H., Shabala S. S., 2013, ApJ, 767, 12
  * Hardcastle (2018) Hardcastle M. J., 2018, MNRAS, 475, 2768
  * Hardcastle & Krause (2013) Hardcastle M. J., Krause M. G. H., 2013, MNRAS, 430, 174
  * Hardcastle & Krause (2014) Hardcastle M. J., Krause M. G. H., 2014, MNRAS, 443, 1482
  * Hardcastle et al. (2019) Hardcastle M. J., et al., 2019, A&A, 622, A12
  * Hurley-Walker et al. (2015) Hurley-Walker N., et al., 2015, MNRAS, 447, 2468
  * Hurley-Walker et al. (2017) Hurley-Walker N., et al., 2017, MNRAS, 464, 1146
  * Hurley-Walker et al. (2019) Hurley-Walker N., et al., 2019, Publ. Astron. Soc. Australia, 36, e048
  * Ineson et al. (2017) Ineson J., Croston J. H., Hardcastle M. J., Mingo B., 2017, MNRAS, 467, 1586
  * Intema et al. (2017) Intema H. T., Jagannathan P., Mooley K. P., Frail D. A., 2017, A&A, 598, A78
  * Joye & Mandel (2003) Joye W. A., Mandel E., 2003, in Payne H. E., Jedrzejewski R. I., Hook R. N., eds, Astronomical Society of the Pacific Conference Series Vol. 295, Astronomical Data Analysis Software and Systems XII. p. 489
  * Jurlin et al. (2020) Jurlin N., et al., 2020, A&A, 638, A34
  * Jurlin et al. (2021) Jurlin N., Brienza M., Morganti R., Wadadekar Y., Ishwara-Chandra C. H., Maddox N., Mahatma V., 2021, A&A, 653, A110
  * Kaiser & Best (2007) Kaiser C. R., Best P. N., 2007, MNRAS, 381, 1548
  * Kakkad et al. (2018) Kakkad D., et al., 2018, A&A, 618, A6
  * Kaviraj et al. (2015) Kaviraj S., Shabala S. S., Deller A. T., Middelberg E., 2015, MNRAS, 452, 774
  * Lacy et al. (2019) Lacy M., et al., 2019, arXiv e-prints, p. arXiv:1907.01981
  * Lynden-Bell (1969) Lynden-Bell D., 1969, Nature, 223, 690
  * Mahatma (2023) Mahatma V. H., 2023, Galaxies, 11, 74
  * Mahatma et al. (2018) Mahatma V. H., et al., 2018, MNRAS, 475, 4557
  * Mahatma et al. (2020) Mahatma V. H., Hardcastle M. J., Croston J. H., Harwood J., Ineson J., Moldon J., 2020, MNRAS, 491, 5015
  * Marshall et al. (2018) Marshall M. A., Shabala S. S., Krause M. G. H., Pimbblet K. A., Croton D. J., Owers M. S., 2018, MNRAS, 474, 3615
  * McConnell et al. (2020) McConnell D., et al., 2020, Publ. Astron. Soc. Australia, 37, e048
  * McNamara & Nulsen (2007) McNamara B. R., Nulsen P. E. J., 2007, ARA&A, 45, 117
  * Mittal et al. (2009) Mittal R., Hudson D. S., Reiprich T. H., Clarke T., 2009, A&A, 501, 835
  * Morganti et al. (2013) Morganti R., Fogasy J., Paragi Z., Oosterloo T., Orienti M., 2013, Science, 341, 1082
  * Morganti et al. (2015) Morganti R., Oosterloo T., Oonk J. B. R., Frieswijk W., Tadhunter C., 2015, A&A, 580, A1
  * Morganti et al. (2021) Morganti R., Oosterloo T., Tadhunter C., Bernhard E. P., Raymond Oonk J. B., 2021, A&A, 656, A55
  * Mukherjee et al. (2021) Mukherjee D., Bicknell G. V., Wagner A. Y., 2021, Astronomische Nachrichten, 342, 1140
  * Murgia et al. (2011) Murgia M., et al., 2011, A&A, 526, A148
  * Myers & Spangler (1985) Myers S. T., Spangler S. R., 1985, ApJ, 291, 52
  * Nesvadba et al. (2008) Nesvadba N. P. H., Lehnert M. D., De Breuck C., Gilbert A. M., van Breugel W., 2008, A&A, 491, 407
  * Norris (2011) Norris R. P., 2011, Journal of Astrophysics and Astronomy, 32, 599
  * Novak et al. (2011) Novak G. S., Ostriker J. P., Ciotti L., 2011, ApJ, 737, 26
  * Parma et al. (2007) Parma P., Murgia M., de Ruiter H. R., Fanti R., Mack K. H., Govoni F., 2007, A&A, 470, 875
  * Pimbblet et al. (2013) Pimbblet K. A., Shabala S. S., Haines C. P., Fraser-McKelvie A., Floyd D. J. E., 2013, MNRAS, 429, 1827
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Predehl et al. (2021) Predehl P., et al., 2021, A&A, 647, A1
  * Quici et al. (2021) Quici B., et al., 2021, Publ. Astron. Soc. Australia, 38, e008
  * Quici et al. (2022) Quici B., Turner R. J., Seymour N., Hurley-Walker N., Shabala S. S., Ishwara-Chandra C. H., 2022, MNRAS, 514, 3466
  * Raouf et al. (2017) Raouf M., Shabala S. S., Croton D. J., Khosroshahi H. G., Bernyk M., 2017, MNRAS, 471, 658
  * Sabater et al. (2019) Sabater J., et al., 2019, A&A, 622, A17
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Shabala & Alexander (2009) Shabala S., Alexander P., 2009, ApJ, 699, 525
  * Shabala & Godfrey (2013) Shabala S. S., Godfrey L. E. H., 2013, ApJ, 769, 129
  * Shabala et al. (2020) Shabala S. S., Jurlin N., Morganti R., Brienza M., Hardcastle M. J., Godfrey L. E. H., Krause M. G. H., Turner R. J., 2020, MNRAS, 496, 1706
  * Shimwell et al. (2019) Shimwell T. W., et al., 2019, A&A, 622, A1
  * Shimwell et al. (2022) Shimwell T. W., et al., 2022, A&A, 659, A1
  * Smith et al. (2016) Smith D. J. B., et al., 2016, in Reylé C., Richard J., Cambrésy L., Deleuil M., Pécontal E., Tresse L., Vauglin I., eds, SF2A-2016: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. pp 271–280 (arXiv:1611.02706)
  * Taylor et al. (2011) Taylor E. N., et al., 2011, MNRAS, 418, 1587
  * Tortora et al. (2009) Tortora C., Antonuccio-Delogu V., Kaviraj S., Silk J., Romeo A. D., Becciani U., 2009, MNRAS, 396, 61
  * Turner (2017) Turner R. J., 2017, PhD thesis, University of Tasmania, Hobart, Australia
  * Turner (2018) Turner R. J., 2018, MNRAS, 476, 2522
  * Turner & Shabala (2015) Turner R. J., Shabala S. S., 2015, ApJ, 806, 59
  * Turner & Shabala (2020) Turner R. J., Shabala S. S., 2020, MNRAS, 493, 5181
  * Turner & Shabala (2023) Turner R. J., Shabala S. S., 2023, Galaxies, 11, 87
  * Turner et al. (2018a) Turner R. J., Rogers J. G., Shabala S. S., Krause M. G. H., 2018a, MNRAS, 473, 4179
  * Turner et al. (2018b) Turner R. J., Shabala S. S., Krause M. G. H., 2018b, MNRAS, 474, 3361
  * Turner et al. (2020) Turner R. J., Drouart G., Seymour N., Shabala S. S., 2020, MNRAS, 499, 3660
  * Turner et al. (2023) Turner R. J., Yates-Jones P. M., Shabala S. S., Quici B., Stewart G. S. C., 2023, MNRAS, 518, 945
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, MNRAS, 444, 1518
  * White et al. (2020a) White S. V., et al., 2020a, Publ. Astron. Soc. Australia, 37, e017
  * White et al. (2020b) White S. V., et al., 2020b, Publ. Astron. Soc. Australia, 37, e018
  * Willott et al. (1999) Willott C. J., Rawlings S., Blundell K. M., Lacy M., 1999, MNRAS, 309, 1017
  * Wolf et al. (2018) Wolf C., et al., 2018, Publ. Astron. Soc. Australia, 35, e010
  * Wright et al. (2010) Wright E. L., et al., 2010, AJ, 140, 1868
  * Yates-Jones et al. (2021) Yates-Jones P. M., Shabala S. S., Krause M. G. H., 2021, MNRAS, 508, 5239
  * Yates et al. (2018) Yates P. M., Shabala S. S., Krause M. G. H., 2018, MNRAS, 480, 5286