# Selecting and Modelling Remnant AGNs with Limited Spectral Coverage

Benjamin Quici1, Ross J. Turner2, Nicholas Seymour1, Natasha Hurley-Walker1
Email<EMAIL_ADDRESS>Stanislav S. Shabala2,3, C. H. Ishwara-
Chandra4
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
3ARC Centre of Excellence for All-Sky Astrophysics in 3 Dimensions (ASTRO 3D)
4National Centre for Radio Astrophysics, TIFR, Post Bag No. 3, Ganeshkhind
Post, 411007 Pune, India

(Accepted XXX. Received YYY; in original form ZZZ.)

###### Abstract

Quantifying the energetics and lifetimes of remnant radio-loud active galactic
nuclei (AGNs) is much more challenging than for active sources due to the
added complexity of accurately determining the time since the central black
hole switched off. Independent spectral modelling of remnant lobes enables the
derivation of the remnant ratio, $R_{\mathrm{rem}}$, (i.e. ‘off-time/source
age’), however, the requirement of high-frequency ($\gtrsim 5$ GHz) coverage
makes the application of this technique over large-area radio surveys
difficult. In this work we propose a new method, which relies on the observed
brightness of backflow of Fanaroff-Riley type II lobes, combined with the
_Radio AGN in Semi-Analytic Environments_ (RAiSE) code, to measure the
duration of the remnant phase. Sensitive radio observations of the remnant
radio galaxy J2253-34 are obtained to provide a robust comparison of this
technique with the canonical spectral analysis and modelling methods. We find
that the remnant lifetimes modelled by each method are consistent; spectral
modelling yields $R_{\mathrm{rem}}=0.23\pm 0.02$, compared to
$R_{\mathrm{rem}}=0.26\pm 0.02$ from our new method. We examine the viability
of applying our proposed technique to low-frequency radio surveys using mock
radio source populations, and examine whether the technique is sensitive to
any intrinsic properties of radio AGNs. Our results show that the technique
can be used to robustly classify active and remnant populations, with the most
confident predictions for the remnant ratio, and thus off-time, in the
longest-lived radio sources ($>50$ Myr) and those at higher redshifts
($z>0.1$).

###### keywords:

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

††pubyear: 2022

##  1 INTRODUCTION

The supermassive black holes (SMBHs) residing at the heart of most galaxies
(Magorrian et al., 1998; Häring & Rix, 2004; Gültekin et al., 2009), play a
profound role in the evolution of their host galaxies and large-scale
(intergalactic) environments. Accretion onto the SMBH powers an active
galactic nucleus (AGN), which, when radio-loud, drives a pair of jets
comprising relativistic plasma. The jets inflate synchrotron-emitting lobes in
the surrounding atmosphere, which act to heat up and expel the surrounding
intergalactic gas, and eventually halt accretion onto the SMBH (e.g. see
reviews by; McNamara & Nulsen, 2007; Fabian, 2012; Alexander & Hickox, 2012).
This _jet mode_ of AGN feedback is needed to quench star formation in the most
massive galaxies (Croton et al., 2006), and to suppress cooling flows in the
cores of massive clusters (Fabian et al., 2003) at low redshifts ($z\lesssim
2-3$; Fabian, 2012).

Prescriptions for these feedback mechanisms require more than just the
energetics and lifetimes of the jetted outbursts (Binney et al., 2007); the
manner by which the jet power couples with the surrounding gas is equally as
important (Tadhunter, 2008; Alexander & Hickox, 2012). The lifetime of the
active jet phase has a large influence on the mechanism coupling the jet power
to the environment; e.g. through shock-heating driven by the global expansion
of the lobes (Worrall et al., 2012), versus the clearing of gas through the
late buoyant rise of jet-inflated bubbles (Churazov et al., 2001). This active
lifetime also strongly maps to the radio source linear size (Kaiser &
Alexander, 1997), which determines the cluster radius of interaction where the
energy deposition takes place. The duration of a quiescent phase, in which
there is either no or reduced jet activity, can additionally be constrained
from the active lifetime when combined with measurements of the duty cycle.
The duty cycle can be directly constrained for “double-double radio galaxies”
(DDRG; Schoenmakers et al., 2000) and some remnant radio galaxies (Turner,
2018), or estimated on a population level from the radio-loud fraction (e.g.
Best et al., 2005; Shimwell et al., 2019). Feedback mechanisms can therefore
be examined by considering the spatial and temporal scales over which the
energy deposition takes place; i.e. does the cluster environment have
sufficient time to reach its equilibrium state before jet activity resumes
(e.g. Turner & Shabala, 2015).

The feedback mechanisms that regulate the supply of the surrounding gas are
thought to give rise to an intermittency in jet activity, and the observed
“remnant” and “restarted” radio galaxy populations (e.g. see reviews by;
Kapinska et al., 2015; Morganti, 2017); the latter representing those in which
jet activity has renewed prior to the remnant lobes fading from view.
Importantly, remnant radio galaxies offer valuable constraints on the active
lifetimes of AGN jets, necessary towards our general understanding of AGN jet
feedback. Shabala et al. (2020) demonstrate that models describing the global
distribution of jet lifetimes and duty-cycles can be inferred from the
observed fractions of remnant and restarted radio sources (relative to active
ones). A direct measurement of the active jet lifetime can be made by
modelling the dynamical evolution of remnant lobes (Turner, 2018); the same is
not true of the age derived for active radio lobes, which corresponds only to
the age of the source at the time of observation, rather than the total
lifetime of the jet. For these reasons, methods to select and model remnant
radio galaxy lobes need to be robust.

Methodologies for selecting remnant radio galaxies are explored by many
authors, particularly since the advent of the Low Frequency Array (LOFAR; van
Haarlem et al., 2013). The overall consensus, as first proposed by Brienza et
al. (2016), is that a combination of morphological- and spectral-based methods
are needed to select a representative sample of remnant radio galaxies.

It is possible to select remnants based on an ultra-steep (e.g. $\alpha\geq
1.2$; Cordey, 1987), or highly curved (e.g. $\Delta\alpha\geq 0.5$; Murgia et
al., 2011) radio spectrum criterion. These methods rely on sufficient fading
(ageing) at the highest observing frequency, and will therefore miss recently
switched-off remnants in which the break frequency lies beyond this frequency;
this is confirmed based on observed and mock radio source populations (e.g
see; Godfrey et al., 2017; Brienza et al., 2017; Mahatma et al., 2018; Shabala
et al., 2020). The bias associated with these spectral-based methods can be
mitigated through the inclusion of a morphological-based selection, where
remnant lobes are selected based on diffuse, low surface-brightness, amorphous
lobes (e.g. Blob 1; Brienza et al., 2016), however such methods may be less
effective for recently switched-off remnants (e.g. 3C 28; Harwood et al.,
2015).

Upon cessation of the jets, a bright radio core is expected to disappear
rapidly ($<0.1\,$Myr). The selection of remnant radio galaxies can therefore
also be made based on the absence of compact radio emission from the core
(e.g; Mahatma et al., 2018; Quici et al., 2021), however we note the recent
results of Jurlin et al. (2021) who demonstrate the co-existence of faint
radio emission from the core together with ultra-steep remnant radio lobes for
a handful of objects. Their results suggest radio emission from the core,
alone, may not necessarily determine whether the radio galaxy is active, and
that high-frequency ($\nu\gtrsim 5$ GHz) observations are needed to confirm
whether radio galaxy lobes are truly remnant.

High frequency radio observations are also critical to modelling the
energetics of both active and remnant AGNs (Turner et al., 2018b). Ageing of
the lobe plasma results in a spectral steepening due to the preferential
radiating of high-energy electrons (Kardashev, 1962; Pacholczyk, 1970), and is
characterised by an optically-thin break frequency that can be related (given
a lobe magnetic field strength) to the plasma spectral age. However, a major
challenge in modelling the dynamical evolution of remnant AGNs comes from the
unknown time since the jets switched off. Turner (2018) show that this problem
can be resolved by modelling the fractional duration of the remnant phase,
based on the observed steepening in the integrated lobe spectrum, reducing the
required model complexity to that of active radio sources. However, this
method requires a well-sampled radio spectrum, where high-frequency
($\nu\gtrsim 5$ GHz) observations are often needed to properly characterise
the observed spectral curvature.

Sensitive, wide-area radio surveys conducted with telescopes such as LOFAR;
the Karoo Array Telescope (MeerKAT; Jonas & MeerKAT Team, 2016); the
Australian Square Kilometre Array Pathfinder (ASKAP; Johnston et al., 2007);
and eventually the Square Kilometre Array (SKA), are expected to offer an
unprecedented view of the radio galaxy life-cycle. In principle, such surveys
enable the compilation of large samples of active, remnant and restarted radio
sources, which will offer invaluable constraints on mechanisms describing the
jet-triggering and AGN feedback. However, it is unrealistic to expect
matching, high-frequency ($\nu\gtrsim 5$ GHz) sky coverage which, at present,
is needed to robustly confirm the status of remnant lobes (Jurlin et al.,
2021), and tightly constrain their energetics (Turner, 2018).

In this work, we therefore explore whether measurements of the off-times of
remnant lobes can be made based on their surface brightness distributions. To
do this, we collect new data, presented in Sect. 2, for the remnant radio
galaxy J2253-34 (see; Quici et al., 2021), which are needed to conduct various
spectral and dynamical modelling of the source. In Sect. 3, we apply an
established spectral ageing method to provide an independent measurement of
the off-time in J2253-34. In Sect. 4, we develop a new dynamical model-based
method that uses the observed lobe backflow to measure the off-times in
remnant radio galaxies. Here, we apply our method to J2253-34, and validate
the estimate of its off-time against that derived previously from its radio
spectrum. In Sect. 7, we use mock radio source populations, simulated at 150
MHz and 1.4 GHz, to examine whether the off-times in a wide range of remnants
can be constrained by this method. Final discussions and conclusions are
summarised in Sect. 6.

The $\Lambda$CDM concordance cosmology with $\Omega_{M}$ = 0.3,
$\Omega_{\Lambda}=0.7$ and $H_{0}$ = 70 km s-1 Mpc-1 (Komatsu et al., 2011) is
assumed throughout the paper.

##  2 DATA

In this section, we describe the observations taken as part of this work, or
larger observing projects, for J2253-34, and the methods used to image these
data. We make use of new data collected by the Australia Telescope Compact
Array (ATCA; Frater et al., 1992), MeerKAT, and the upgraded Giant Meterwave
Radio Telescope (uGMRT), together with several ancillary data products. A
multi-wavelength view of J2253-34 using these new data can be seen in Fig. 1.
All measurements of the integrated flux density are reported in Tab. 9.

Figure 1: The remnant radio galaxy J2253-34 and its intergalactic environment.
The background image is an RGB composite constructed from mid-infrared $K_{s}$
(red), $H$ (green) and $J$ (blue) bands, respectively (e.g. see; Quici et al.,
2021). The host galaxy is indicated by a cyan-coloured ‘x’ marker. A 5.5 GHz
view of J2253-34 with $4^{\prime\prime}$ resolution is superimposed onto the
background. Coloured contours show the 1.42 GHz radio emission (‘sub-band 3’
in Sect. 2.2), with levels increasing uniformly in log-space by $0.2$ dex
between [$3,80$] $\times\sigma$. Ellipses in the lower-left corner denote the
shape of the restoring beam at 1.42 GHz (orange) and 5.5 GHz (blue).

###  2.1 ATCA data

Table 1: Summary of the new ATCA observations in different configurations described in Sect. 1. In order, columns denote: the principal investigator; the project code; the ATCA configuration; the date of observation, each conducted during the year 2021; the total integration time; the range of shortest to longest baseline lengths; and the largest recoverable scales at 5.5 GHz and 9 GHz. We note, Bmax quoted for all but the 6C configuration excludes baselines formed with the $\sim$6 km outrigger antenna (CA06), due to the gap in $(u,v)$ coverage. PI | Project code | Config. | Date | $t$ | Bmin – Bmax | $\theta_{\mathrm{max,5.5}}$ | $\theta_{\mathrm{max,9}}$
---|---|---|---|---|---|---|---
dd/mm | (h) | (m) | (′) | (′)
Benjamin Quici | C3335 | H75 | $03/06$ | 8 | 31 – 89 | 6.0 | 3.7
C3335 | H168 | $22/07$ | 5.5 | 61 – 185 | 3.1 | 1.9
C3402 | 750C | $14/02$ | 12 | 46 – 704 | 4.1 | 2.5
C3402 | 1.5A | $03/01$ | 12 | 153 – 1469 | 2.4 | 1.5
C3402 | 6C | $(12,15)/11$ | 10 | 337 – 5939 | 1.2 | 0.7

Pointed observations of J2253-34 were awarded through a series of observing
proposals (Tab. 1). Observations were conducted within the $4$ cm–band, which
allows simultaneous observing over two 2 GHz-wide bands which were centered at
5.5 GHz and 9 GHz. The collection and reduction of these data are presented
below.

To produce images with high spatial-resolution and sensitivity to large-scale
structure, coverage of the $(u,v)$ plane was maximised by observing in various
ATCA configurations (see Tab. 1). Baselines ranged between 31–5939 m, ensuring
an appropriate sampling of the largest angular scale of J2253-34 ($\theta\sim
100^{\prime\prime}$). The combined integration time on the target was 39.5
hours. At the beginning of each observation, we observed PKS B1934-638 for 10
minutes to derive bandpass, gain and flux calibration solutions. 12 minute
scans of J2253-34 were interleaved with two-minute scans of PKS B2255-282; a
source used to calibrate the time-varying gains.

Data reduction was performed using the Multichannel Image Reconstruction,
Image Analysis and Display (miriad; Sault et al., 1995), and closely followed
the method outlined by Huynh et al. (2015, 2020); we used the primary scan to
solve for the shape of the bandpass, the secondary scan to calibrate for the
time-varying gains and the polarization leakages, and then transferred the
solutions to J2253-34 where the gain solutions were averaged over a two-minute
interval (the length of each phase calibrator scan). Throughout all these
steps, channels known to be contaminated by radio frequency interference (RFI)
were automatically flagged.

Calibrated data-sets associated with each unique observation, with equal
weighting, into a single data-set. To improve the imaging dynamic range, three
rounds of phase-only self calibration were conducted using a Robust 0.0 Briggs
weighting (Briggs, 1995), corresponding to an image resolution of
approximately $3^{\prime\prime}\times 6^{\prime\prime}$. For each round, gains
were solved over 2, 1, and 0.5 minute time intervals, respectively. We refer
to this data below as the ‘calibrated data-set’, and discuss the synthesis of
our science-ready images in Sects. 2.1.1 & 2.1.2. For these data, we assume a
$3\%$ calibration uncertainty in the ATCA flux density scale (Reynolds, 1994).

####  2.1.1 Low-resolution imaging

We first produced images with lower spatial resolution to improve constraints
on the high-frequency, integrated, spectral energy distribution (SED). This
was done by splitting each 2-GHz band of the calibrated data-set into four 512
MHz sub-bands. Each sub-band was imaged in miriad using a Robust 0.0 Briggs
weighting, with an additional $30^{\prime\prime}$ Gaussian taper on the
visibilities so as to reduce contribution and noise from the longer baselines.
The data were CLEANed to an initial depth of 3$\sigma$, after which a deeper
CLEAN to 1$\sigma$ was performed around pixels already in the model. Images
were corrected for attenuation of the primary beam.

To extract the photometry, we used
polygon_flux111https://github.com/nhurleywalker/polygon-flux, which uses the
shape of the restoring beam to evaluate a background-subtracted integrated
flux density corresponding to a user-defined region (e.g. see Hurley-Walker et
al. 2019). Measurement uncertainties were calculated by multiplying the RMS
noise by the square-root number of beams corresponding to the region volume,
which we added in quadrature with the absolute calibration uncertainty.

####  2.1.2 High-resolution imaging

We further produced images with high spatial resolution so we could map the
radio spectrum across the lobes (in conjunction with our data at other
frequencies). To maximise signal-to-noise, we imaged the entire 2 GHz
bandwidth available in each wide-band. To produce these images, we used a
uniform Briggs weighting and applied a $8^{\prime\prime}$ taper to the
visibilities; although it was possible to improve the spatial resolution
beyond $\theta_{\mathrm{FWHM}}=8^{\prime\prime}$, this corresponded to the
best possible circular resolution that could be achieved across all radio
images presented in this work. To demonstrate this in Fig. 1, we present a
higher-resolution view of J2253-34 at 5.5 GHz by instead applying a
$4^{\prime\prime}$ taper to the visibilities. The data were CLEANed to an
initial depth of 3$\sigma$, after which a deeper CLEAN to 1$\sigma$ was
performed around pixels already in the model. After correcting images for the
frequency-dependent primary beam attenuation, they achieved a mean root mean
square (RMS) sensitivity of 11 $\mu$Jy beam-1 and 13 $\mu$Jy beam-1 at 5.5 GHz
and 9 GHz, respectively. We refer to these images as ‘high-resolution ATCA
images’.

###  2.2 MeerKAT data

We obtained 3.5 hours of time to observe J2253-34 on the MeerKAT radio
telescope (PI: B. Quici, project code: MKT 20212). Observations took place in
the ‘L-band’ (0.89–1.71 GHz), thus offering sensitive, high spatial-resolution
observations at an intermediate frequency range. To collect the data, we used
PKS 1934-638 to calibrate the bandpass, and PKS 2259-37 to calibrate the time-
varying gains. The MeerKAT bandpass response is stable to within 3% over 2–3
hours. As such, the observations were conducted in two ‘blocks’ by visiting
the bandpass calibrator twice; once at the beginning of the observation, and
again 1.5 hours later. Each bandpass calibrator scan lasted 10 minutes.
Observations of J2253-34 were carried out by interleaving a two minute gain
calibrator scan with a 15.5 minute target scan. In total, the target received
2.6 hours of exposure time. Using the calibrated measurement set made
available from the South African Radio Astronomy Observatory (SARAO) Web
Archive222The SARAO Web Archive is found at https://archive.sarao.ac.za/, we
performed a single round of RFI flagging using AOFlagger (Offringa et al.,
2012). We identified four frequency ranges as having RFI-clean, unflagged
visibilities, and we refer to these as ‘sub–bands’ (see Tab. 9).

Imaging was conducted using WSClean (Offringa et al., 2014). Each sub–band was
imaged independently by blindly CLEANing the data down to a depth of
3$\sigma$, followed by a deeper CLEAN down to a $1\sigma$ level around the
model. Imaging of sub–bands 1–4 was performed using an image weighting of
robust -1.0, -0.6, 0.0 and +0.25, respectively. In this way, all but the
lowest sub–band achieved an image resolution of
$\theta_{\mathrm{FWHM}}\lesssim 8^{\prime\prime}$. All imaging was performed
using the multiscale option, which enables multi-scale deconvolution (Offringa
& Smirnov, 2017). Each sub–band was then corrected for the frequency-dependent
primary beam attenuation. Finally, images associated with sub–bands 2–4 were
convolved to a circular 8′′ resolution, in order to match the resolution in
the high-resolution ATCA images. The image associated with sub–band 1 was
excluded from this final step, as the restoring beam exceeded
$8^{\prime\prime}$. A 5% calibration uncertainty is assumed throughout this
work.Following the method as described in Sect. 2.1.1, the integrated flux
density of J2253-34 (and associated measurement uncertainty) was measured from
each sub–band using the same user-defined polygon.

###  2.3 uGMRT data

We finally made use of two low-frequency, high-resolution radio images,
created using data collected with the upgraded Giant Meterwave Radio Telescope
(uGMRT), which contained J2253-34 within their fields of view. Observations
were conducted in band-3 ($250-500$ MHz; proposal code 37_057, PI: Ross
Turner), and in band-4 ($550-900$ MHz; proposal code 35_022, PI: C. H. Ishwara
Chandra). The data was reduced using the CAPTURE Pipeline (Kale & Ishwara-
Chandra, 2021), which uses the Common Astronomy Software Applications (casa;
McMullin et al., 2007) package. After the initial gain calibration, the data
was channel averaged to keep the bandwidth smearing negligible. Using a Briggs
weighting of robust 0.0, the channel averaged source data was subjected to
four rounds of phase-only self-calibration followed by 4 more rounds of
amplitude and phase self-calibration. A final band-3 image, centered at $417$
MHz, was produced using a robust 0.0 weighting, resulting in a native
resolution of $\sim 8^{\prime\prime}\times 5^{\prime\prime}$. A final band-4
image, centered at $682$ MHz, was produced using a robust 0.0 weighting,
resulting in a native resolution of $\sim 5^{\prime\prime}\times
3^{\prime\prime}$. The resolution in each image was convolved to a circular
beam of $\theta=8^{\prime\prime}$, achieving a local RMS of $52\mu$Jy beam-1
and $43\mu$Jy beam-1 at 417 MHz and 682 MHz, respectively.

###  2.4 Ancillary data

To fill out the sampled integrated radio spectrum of J2253-34, particularly at
frequencies not covered by our new data, we take several measurements
previously compiled by Quici et al. (2021). We take their integrated flux
densities reported at: 119 and 154 MHz, measured from the GaLactic and
Extragalactic All-sky Murchison Widefield Array (GLEAM; Hurley-Walker et al.,
2017) survey; 887 MHz, measured from the Evolutionary Map of the Universe
(EMU; Norris, 2011); 1.4 GHz, measured from the NRAO VLA Sky Survey (NVSS;
Condon et al., 1998), and; 2 and 2.868 GHz, measured using pointed
observations collected in the $LS$-band ($\lambda=4$ cm) of ATCA (see Tab. 9).

##  3 SPECTRAL MODELLING OF REMNANT LOBES

Before we develop a new method to measure the off-time in remnant radio
sources, we must ensure that our result can be verified against an established
method. The off-time can be derived from the rapid steepening of their high-
frequency spectrum due to the cessation of the jet activity and the subsequent
lack of freshly shock-accelerated synchrotron-emitting electrons. In Sect.
3.1–3.3, we outline the theory of several literature-established spectral
ageing models, and discuss their applications to radio lobe spectra. In Sect.
3.4, these models are applied to a combination of low- and high resolution
radio observations, in order to tightly parameterise the spectral ageing
observed for J2253-34. Further, in Sect. 3.2, we briefly describe a python
package we developed as part of this work to model synchrotron spectra for a
diverse range of applications in astrophysics.

###  3.1 Synchrotron emissivity

The total emissivity of a synchrotron-emitting plasma, $J(\nu)$, can be found
by integrating the single-electron emissivity (e.g. $j(\nu)$; Longair, 2011)
over the distribution of electron energies, $N(E)$; and the probabilistic
distributions of the magnetic field strength, $p_{B}$, and pitch angle,
$p_{\xi}$ (Equation 4 of; Hardcastle, 2013)

$J(\nu)=\int_{0}^{\infty}\int_{0}^{\pi}\int_{0}^{\infty}\frac{\sqrt{3}Be^{2}\sin\xi}{8\pi^{2}\epsilon_{0}cm_{e}}F(x)N(E)p_{\xi}p_{B}dEd\xi
dB$,

(1)

where $B$ gives the magnetic field strength; $\xi$ the pitch angle; $c$ the
speed of light; $e$ the electron charge; $\epsilon_{0}$ the permittivity of
free space; and $m_{e}$ the electron mass. Rybicki & Lightman (1979) define
$F(x)$ as the shape of the single-electron radiation spectrum, where $x$
represents a dimensionless parameterisation of frequency, magnetic field
strength, and electron energy, $E$, through:

$x=\frac{\nu}{\nu_{\mathrm{b}}}=\frac{4\pi
m_{\mathrm{e}}^{3}c^{4}\nu}{3eE^{2}B\sin\xi},$ (2)

where $\nu_{\mathrm{b}}$ is the break frequency above which (i.e. at higher
frequencies) the spectrum steepens due to energy losses associated with the
plasma ageing. Importantly, for a given magnetic field strength, the break
frequency is related to the plasma spectral age through:

$\tau=\frac{\upsilon
B^{1/2}}{B^{2}+B_{\mathrm{IC}}^{2}}\big{(}\nu_{\mathrm{b}}(1+z)\big{)}^{-1/2},$
(3)

where $B_{\mathrm{IC}}=0.318(1+z)^{2}$ nT gives the magnetic field equivalent
of the cosmic microwave background at redshift $z$, and $\upsilon$ is a
constant of proportionality (see Equation 4 of Turner et al. 2018b).

Parameterisation of meaningful attributes, such as $\nu_{\mathrm{b}}$, demands
appropriate treatment of $N(E)$, $p_{B}$, and $p_{\xi}$. Kardashev (1962)
initially proposed a model with no pitch angle scattering, giving $p_{\xi}=1$,
however a more realistic model assuming an isotropic distribution for the
pitch angle was subsequently suggested by Jaffe & Perola (1973), giving
$p_{\xi}=\frac{1}{2}\sin\xi$; we assume the latter probabilistic distribution
for the pitch angle in this work. Meanwhile, it is common to assume a locally-
uniform magnetic field, such that $p_{B}=1$, however, this assumption is
likely not physical due to the turbulent mixing of lobe plasma. Tribble (1991)
provide an alternative treatment for a turbulent (non-uniform) magnetic field,
in which the field strength is drawn from a Maxwell-Boltzmann distribution;
here, $p_{B}$ is satisfied by

$p_{B}=\sqrt{\frac{2}{\pi}}\frac{B^{2}\exp(-B^{2}/2a^{2})}{a^{3}},$ (4)

where $a=B_{0}/\sqrt{3}$ for mean magnetic field strength, $B_{0}$
(Hardcastle, 2013). These assumptions can be paired together with an
appropriate treatment of the electron energy distribution of the underlying
plasma, $N(E)$ (described below), to create a range of spectral ageing models;
these are referred to as either the “standard” (uniform $B$) or “Tribble”
(Maxwell-Boltzmann distribution $B$) forms of each spectral ageing model.

####  3.1.1 Impulsive injection models

The radio spectrum arising from an impulsively-injected electron population is
modelled as an ensemble of equal-age electrons [$t,t+dt$] injected at $t=0$
(Kardashev, 1962; Jaffe & Perola, 1973). Electrons are injected with an
initial power-law distribution of energies, e.g $N(E)=N_{0}E^{-s}$, where the
power-law energy index, $s$, is related to the injection spectral index,
$\alpha_{\mathrm{inj}}$, through $s=2\alpha_{\mathrm{inj}}+1$. Synchrotron and
inverse-Compton (IC) radiative losses, the latter due to up-scattering of
Cosmic Microwave (CMB) photons, occur as the electron packet ages; the
resulting effect on the energy distribution is well understood (e.g. see;
Pacholczyk, 1970; Longair, 2011). In the case of a uniform magnetic field, we
follow the method of Turner et al. (2018b) to remove the dependence on $B$ by
recasting $N(E)$ to:

$N(x)=x^{-1/2}\begin{cases}(x^{1/2}-\iota^{1/2})^{s-2}&x>\iota\\\
0&x<\iota\end{cases}$ (5)

where $\iota(\nu,\xi)=\nu/(\nu_{\mathrm{b}}\sin\xi)$ assuming an isotropic
distribution for the pitch angle, as proposed by Jaffe & Perola (1973). Here,
the total synchrotron emissivity is thus evaluated as:

$\mathcal{J}(\nu)=J_{0}\,\nu^{(1-s)/2}\int_{0}^{\pi/2}\sin^{(s+3)/2}\xi\int_{0}^{\infty}F(x)N(x)~{}dx\,d\xi,$
(6)

where $J_{0}$ acts as a frequency-independent constant. In this work, we refer
to this impulsively-injected spectral ageing model as the “standard” Jaffe-
Parola (or JP) model. Importantly, parameterising the shape of this standard
JP spectrum in this way only requires knowledge of the injection index and
break frequency. By comparison, the shape of a “Tribble” Jaffe-Parola (or TJP)
model spectrum additionally requires knowledge of $B$.

Considering some small region across a radio galaxy lobe (e.g. bounded by the
FWHM of a restoring beam), it is reasonable to suspect that the plasma
contained here will have been injected at approximately the same age. This
means that the associated radio spectrum is likely that of the JP model,
however we note the results of Turner et al. (2018a), Mahatma et al. (2020)
and Yates-Jones et al. (2021, 2022) who show that mixing of plasma ages may
violate this assumption.

####  3.1.2 Continuous Injection models

The above assumption of an impulsively-injected packet of electrons breaks
down when considering the integrated spectrum of a radio AGN, which comprises
multiple, or on-going, such injections throughout the lifetime of the source.
A continuous injection model describes a radio spectrum arising from a time-
averaged population of continuously-injected electron packets; that is, a
mixed-age plasma population with ages ranging within [$0,\tau$]. In the
simpler and more common case, the CI model assumes continuous injection via an
active jet throughout the entire lifetime of the source (e.g, the ‘CI-on’ or
CIJP333Here, the ‘JP’ reflects only the isotropic distribution for the pitch
angle proposed by Jaffe & Perola (1973), not their assumption of an
impulsively-injected plasma. model). To model remnant AGNs, Komissarov &
Gubanov (1994) present an extension to the CI model, known as the ‘CI-off’ or
KGJP model, by instead assuming some fraction of the total source age,
$R_{\mathrm{rem}}=t_{\mathrm{rem}}/\tau$, is spent in a remnant phase; we
refer to $R_{\mathrm{rem}}$ as the remnant ratio, where the total source age,
$\tau$, is given by the sum of the active and remnant durations,
$\tau=t_{\mathrm{on}}+t_{\mathrm{rem}}$. As with the JP model, the effect of
ageing on the plasma energy distribution is well understood (e.g. see;
Longair, 2011). In the case of a uniform magnetic field, the recast form of
$N(E)$ is given by (Equation 3 of; Turner, 2018)

$\mathcal{N}(x)=x^{-1/2}\begin{cases}\leavevmode\resizebox{}{7.2pt}{$(x^{-1/2}-\zeta^{1/2})^{s-1}-(x^{1/2}-\iota^{1/2})^{s-1}$}&x>\iota\\\
\leavevmode\resizebox{}{7.2pt}{$(x^{1/2}-\zeta^{1/2})^{s-1}$}&\zeta\leq
x\leq\iota\\\ 0&x<\zeta\end{cases},$ (7)

where $\zeta(\nu,\xi)=\iota(\nu,\xi)R_{\mathrm{rem}}^{2}$. The total
synchrotron emissivity of a CI spectrum is thus derived by:

$\mathcal{J}(\nu)=J_{0}\,\nu^{-s/2}\int_{0}^{\pi/2}\sin^{(s+4)/2}\xi\int_{0}^{\infty}F(x)N(x)~{}dx\,d\xi,$
(8)

In the limiting case where $R_{\mathrm{rem}}=0$, the CI-off model reduces to
the simpler CI-on form. For simplicity, we thus refer to the general form as
the CI model. In this way, parameterising the shape of a standard CI spectrum
requires only the injection index, break frequency, and remnant ratio. As with
the TJP model, the shape of a “Tribble” Continuous Injection (or TCI) spectrum
additionally requires knowledge of $B$.

Due to their assumption of a mixed-age plasma, the CI models are relevant for
deriving the plasma spectral age from the integrated spectra of radio lobes
(e.g; Parma et al., 2007; Murgia et al., 2011; Shulevski et al., 2012; Brienza
et al., 2016; Turner, 2018; Duchesne & Johnston-Hollitt, 2019). However, it is
worth mentioning that the CI model assumes no temporal evolution in the
magnetic field, which we know not to be true for expanding lobes. Using
synthetic radio source populations, for which all radiative losses and
temporal evolution of the magnetic field were considered, Turner et al.
(2018b) demonstrated that the radio emissivity reflects only the most recent
$\sim 10\%$ of the magnetic field history. Importantly, their results showed
that the spectral age derived using the break frequency, parameterised by the
CI model, together with the present-time magnetic field strength, was in
agreement with the dynamical age. As such, we can be confident that the
quantities parameterised by the CI model are physically meaningful.

###  3.2 Spectral model fitting

We now discuss our approach taken to numerically implement the spectral ageing
model described in the previous section, and to fit these spectra to observed
data points. The code we developed to achieve this is expected to be highly
valuable to the broader astrophysics community (e.g. synchrotron spectra from
active galaxies, supernovae remnants, etc.), and so we have released it as the
publicly-available synchrotron spectral fitter
(synchrofit)444https://github.com/synchrofit python library. Principally,
synchrofit offers a suite of routines that: 1) implement the JP and CI
spectral ageing models described in Sect. 3.1.1 and 3.1.2, and 2) fit these
models to real radio data.

For computational efficiency, the standard spectral ageing models are
evaluated by numerically integrating over the electron energy and pitch angle,
following the reduced dimensionality parameterisation of $J(\nu)$, first
presented by Nagai et al. (2006, JP), Turner et al. (2018b, CI-on) and Turner
(2018, CI-off). The Tribble spectral ageing models are evaluated in a similar
fashion, however involve an additional numerical integration over the magnetic
field strength, and as such incur a considerable cost to the computation time.
An injection index and break frequency are required to parameterise the shape
of the JP spectrum, with an additional normalisation to scale the fit
appropriately. The shape of a CI spectrum is additionally parameterised by the
remnant ratio, noting that $R_{\mathrm{rem}}=0$ in the case of the CI-on
model. Tribble forms of the spectral ageing models are parameterised
similarly, however additionally require the redshift and magnetic field
strength.

Fitting of these models is conducted using a grid-search informed by an
adaptive maximum likelihood algorithm. Discrete samples in $s$,
$\log{\nu_{b}}$ and $R_{\mathrm{rem}}$ are created by uniformly sampling their
allowed ranges, which are used to construct a three-dimensional parameter
grid. A spectral model is evaluated for each unique parameter combination,
where a corresponding normalisation is optimised following
$S_{0}=\sum^{n}_{i=1}(S_{\mathrm{o_{i}}}\times\sigma_{i})-\sum^{n}_{i=1}(S_{\mathrm{p_{i}}}\times\sigma_{i})/\sum^{n}_{i=1}\sigma_{i}$,
where $n$ is the number of spectral measurements, $S_{\mathrm{o_{i}}}$ and
$S_{\mathrm{p_{i}}}$ are the observed and predicted flux densities in
logarithmic space, and $\sigma_{i}$ is the measurement uncertainty for the
$i$th datum. The quality of fit is evaluated using the lowest value of the
Akaike information criterion (AIC; Akaike, 1974) following
$\text{AIC}=2k-2\ln\mathcal{L}$, where $k$ is the number of model
parameters555For the standard JP and CI-on models, $k=3$. For the standard CI-
off models, $k=4$. In the case of the Tribble forms, $k$ increases by 1 due to
the magnetic field.. The likelihood function, $\mathcal{L}$, is defined
through:

$\mathcal{L}(\mu,\sigma)=\prod_{i=1}\frac{1}{\sqrt{2\pi}\sigma_{i}}\exp\Big{(}-\frac{(x_{i}-\mu_{i})^{2}}{2\sigma_{i}^{2}}\Big{)}$
(9)

which, for each unique measurement, $i$, takes into consideration the observed
value, $x_{i}$, measurement uncertainty, $\sigma_{i}$, and the corresponding
predicted model value $\mu_{i}$. Across all grid realisations, the set of
parameters providing the optimal spectral fit are defined as those that
minimise the AIC. The relative quality of each spectral fit is found using
$p_{i}=e^{(\mathrm{AIC}_{i}-\mathrm{AIC}_{\mathrm{best}})/2}$, where $p_{i}$
is the probability that the $i^{th}$ spectral model provides a better fit to
the observations than the selected ‘best’ model. Provided the marginal
probability density functions (PDFs) in each model parameter are pseudo-
Gaussian, the fitted uncertainty in each parameter is taken as the standard
deviation, and the peak in the PDF as the approximate location of the global
maximum. In order to increase the precision of the best estimate in each free
parameter, a grid of points is then re-sampled around the global maximum; in
this work, we repeat this over three iterations. We note, Gaussianity may not
necessarily always hold for these PDFs, as this will depend upon the spectral
fit, hence we make sure to inspect these throughout this work.

###  3.3 Effect of inhomogeneity in the local magnetic field

Figure 2: Comparison between standard (solid lines) and Tribble forms (broken
lines) of the JP (top panel) and CI (bottom panel) spectra. synchrofit is used
to simulate radio spectra for $s=2.2$ and $\nu_{\mathrm{b}}=2.1$ GHz. Tribble
models are shown for $B=0.1$ nT (dot-dashed), $B=1$ nT (dotted), and $B=10$ nT
(dashed). CI models are shown for $R_{\mathrm{rem}}=0$ (orange) and
$R_{\mathrm{rem}}=0.5$ (purple). We note that the Tribble models almost
exactly overlay each other and hence can not be distinguished.

The standard form of the spectral ageing models have a considerable
computational advantage over the (likely more physically accurate) Tribble
forms in our implementation (due to the lower number of numerical
integrations). We therefore investigate if the standard form can be used to
reasonably approximate the Tribble form for either spectral ageing model, or
for any range of magnetic field strengths.

Both the standard and Tribble forms of the spectral ageing models offer a
physically-meaningful parameterisation of observed radio spectra. Due to the
consideration of a variable magnetic field, the shape of the expected
synchrotron spectrum differs to that in which the magnetic field is uniform.
For example, Hardcastle (2013) show a significant departure of the TJP model
from the standard form, notably due to the sharper spectral turnover around
the break frequency in the TJP spectrum. Turner et al. (2018b) explore the
consequence of using the standard spectral ageing models to parameterise
simulated Tribble spectra; their results demonstrate an over-estimation of the
break frequency, which appears to be more appreciable for the JP model as
compared to the CI model.

The assumption of a locally uniform magnetic field will be at least somewhat
violated in radio AGN lobes due to turbulent flow; the Maxwell-Boltzmann
magnetic field distribution of the Tribble spectral ageing models may provide
a more realistic description, however this in turn requires knowledge of
$B_{0}$ which, unless independent X-ray data is available, would require
dynamical modelling.

As such, we seek to explore the dependence of the shape of the Tribble
spectral ageing models on $B_{0}$. We use synchrofit to simulate standard and
Tribble forms of the JP and CI models, assuming $s=2.2$, and
$\nu_{\mathrm{b}}=2.1$ GHz. Tribble models are simulated for $B_{0}=0.1$ nT,
$1$ nT and $10$ nT, in order to cover a broad range of acceptable values.
These are shown in Fig. 2.

As a consistency check with Hardcastle (2013), Harwood et al. (2013) and
Turner et al. (2018b), we reassuringly find that integrating over a Maxwell-
Boltzmann magnetic field modifies the shape of the predicted spectra to that
in which the magnetic field is uniform. However, our results show that the
shape of the predicted Tribble spectra are largely insensitive to values of
the mean magnetic field strength. We find that electron packets with magnetic
strengths in the exponential tail of the Maxwell-Boltzmann distribution (i.e.
$B>\sqrt{{3}/{2}}\,B_{0}$) contribute the vast majority of the total
synchrotron emission at all radio frequencies. Importantly, the mean magnetic
field strength has no effect on the shape of this exponential portion of
probability density function, albeit modifies its normalisation. As a result,
the shape of the radio spectrum remains largely unaffected for moderate
changes in the mean magnetic field.The implication here is it is possible to
accurately parameterise the injection index, break frequency, and remnant
fraction (in the case of the CI models), using the Tribble spectral ageing
models without explicitly requiring knowledge of $B_{0}$; an informed guess of
the mean magnetic field strength is sufficient (e.g. $B_{0}\sim 0.1$ \- $10$
nT).

###  3.4 Modelling the radio spectrum of J2253-34

Figure 3: The injection index optimised for J2253-34. Panels a) and b) show
the observed radio emission at 1.32 GHz and 5.5 GHz, respectively; each image
is presented on a square-root stretch, and is normalised between the
minimum/maximum pixel values. A solid dark-orange contour denotes the
5$\sigma$ brightness of the source at 9 GHz. Circular markers of an 8′′
aperture (the FWHM of the restoring beam) are used to represent independent
regions from which a radio spectrum is extracted and used as a constraint on
the injection index. Apertures overlaid on the south-western lobe are colour-
coded and numbered based on aperture number, and are highlighted in panels c)
and d); to avoid clutter in the figure we do not plot the spectral fits for
apertures in the north-eastern lobe. Panel c) demonstrates fits of the TJP
model ($B_{0}=0.1$ nT) optimised to their associated spectrum. Panel d)
demonstrates the probability density function associated with each spectral
fit, marginalised with respect to the injection index; these are also
represented by the box and whisker plots, where the box is constrained by the
$25^{\mathrm{th}}$ \- $75^{\mathrm{th}}$ percentiles. A global optimisation
for the injection index is found by multiplying the PDFs associated with each
lobe (dot-dashed black curve), which is then approximated by a 1D Gaussian
(solid red curve) in order to characterise the peak probable value and
uncertainty: $s~{}=~{}2.25\,\pm\,0.03$ (or alternatively, the injection
spectral index; $\alpha_{\mathrm{inj}}=0.625\,\pm\,0.015$).

We are now equipped with the tools needed to derive a robust measurement of
the off-time from our observations J2253-34. However, we choose not to
directly fit the integrated spectrum with the CI-off model as the three
parameters are somewhat correlated, and may be misfitted if our observations
include moderate systematic uncertainties. Instead, we independently constrain
the injection index by considering the spectra arising from narrow regions
across the lobes; it is reasonable to assume these spectra can be described by
the JP model (see Sect. 3.1.1). Assuming a constant injection index throughout
the lifetime of the source, the JP models can be jointly fit to these spectra,
resulting in a precise measurement of the injection index optimised across the
lobes. The integrated spectrum can now be fitted using the CI model for a
known injection index, reducing any potential systematic errors in estimates
for the break frequency and remnant ratio; the latter mapping directly to the
duration of the remnant phase for any given source age. This is described in
detail below.

Figure 4: The spectral emissivity of J2253-34. Blue markers are used to
display measurements of the integrated flux density. Error-bars represent
1$\sigma$ measurement uncertainties. Overlaid is the fit of the TCI model,
assuming $B_{0}=0.1\,$nT and $s=2.25\pm 0.03$ (Sect. 3.4.1). The optimal fit
of the TCI model (Sect. 3.4.2) is displayed using a solid black line. Also
shown in light-gray are all plausible spectral fits within a 99.7% ($3\sigma$)
confidence interval. The spectral break frequency and remnant ratio are fitted
as $\nu_{\mathrm{b}}~{}=~{}1.26_{-0.05}^{+0.06}$ GHz, and
$R_{\mathrm{rem}}=0.23\pm 0.02$, respectively.

####  3.4.1 Modelling the injection index

To model the injection index, we use the 8′′ beam-matched radio maps described
in Sect. 2. We use miriad to create a 2D cutout of J2253-34 from each image,
and regrid the cutouts to have a common pixel scale
($d\theta=1^{\prime\prime}$) and projection (SIN). To efficiently sample the
resolved radio spectrum, we extract the radio spectra arising from individual
circular apertures arranged in a hexagonal configuration across the lobes. The
full width half maximum (FWHM) of each aperture is set to
$\theta_{\mathrm{FWHM}}=8^{\prime\prime}$ to reflect the size of the restoring
beam. The integrated flux density within each aperture is then calculated
following $S_{\mathrm{int}}=\sum^{n}_{i=1}S_{p,i}\times d\theta^{2}/\Omega$,
where $n$ is the number of pixels contained within each aperture, $S_{p,i}$ is
the value of the $i^{\mathrm{th}}$ pixel in Jy beam-1, $d\theta$ is the pixel
angular scale, and $\Omega$ is the beam solid angle defined as
$\theta_{\mathrm{FWHM}}^{2}\times\pi/(4\ln 2)$. Given the spectra are
extracted from exactly one resolution element, the measurement uncertainty is
evaluated as the quadrature sum of the RMS noise in the observations and their
absolute calibration uncertainty.

Following this method, the radio spectra are extracted from all regions across
the lobes with 9 GHz radio emission above the 5$\sigma$ noise level. This
restriction is done so as to maximise the number of observing frequencies
contributing to the radio spectrum. The spectra arising from each aperture are
independently fit using the JP (and TJP) models; importantly both the break
frequency and the injection index are free parameters in their corresponding
aperture spectral fit. In light of the results in Sect. 3.3, we use Equation 2
of Miley (1980) to calculate a crude approximation of the lobe equipartition
magnetic field strength as $B_{\mathrm{eq}}\approx 0.35\,$nT. Considering also
that radio galaxy lobes are typically sub-equipartition (e.g. see Fig. 11 of;
Turner et al., 2018b), we fit the TJP model assuming $B_{0}=0.1$ nT, $0.35$ nT
and 1 nT. Spectral models are evaluated by uniformly sampling the injection
index within a range corresponding to [2.01, 2.99] and an initial grid spacing
of $\Delta s~{}=~{}0.05$. For each aperture spectral fit, a corresponding PDF
is found by marginalising across the fitted injection index. Each PDF is then
multiplied together to determine which value for the injection index provides
an optimal fit across all aperture spectra; this assumes no correlation
between the spectra arising from adjacent apertures, and also assumes a
constant injection index throughout the lifetime of the source. Following this
method, we optimise an injection index by considering the spectra arising from
both lobes, which we present in Fig. 3. Importantly, we found that the optimal
injection index modelled using the standard JP model is found to be
$s~{}=~{}2.10\,\pm\,0.05$ (or; $\alpha_{\mathrm{inj}}=0.550\,\pm\,0.025$). For
the Tribble model, we found an optimal value of $s~{}=~{}2.25\,\pm\,0.03$ (or;
$\alpha_{\mathrm{inj}}=0.625\,\pm\,0.015$).

####  3.4.2 Modelling the break frequency & remnant ratio

We compile an integrated radio spectrum for J2253-34 using the two decades of
radio-frequency data described in Sect. 2. The raw spectrum is presented in
Fig. 4, and demonstrates a clear spectral curvature owing to the losses
associated with an inactive jet; this is evidenced by an ultra-steep high-
frequency spectral index of $\alpha=1.7$, computed between 4.5 - 9.5 GHz. With
the injection index estimated in Sect 3.4.1, we seek to model the break
frequency and remnant ratio of J2253-34 by parameterising the observed
curvature using the CI spectral ageing models.

Similar to the resolved modelling to determine the injection index, we model
the integrated spectrum using both the standard and Tribble forms of the CI
model. For consistency, we use the injection indices estimated using the
standard and Tribble JP models as inputs to the standard and Tribble CI
models, respectively. Similarly, TCI models are evaluated for the same values
of $B_{0}$ used in Sect. 3.4.1. To account for the small uncertainty in the
injection index previously fitted by the JP models, we treat ‘$s$’ as a free
parameter in the CI models, but apply a tight Gaussian prior based on the
fitted uncertainty. We find that fixing the injection index offers no
appreciable difference to our results, unsurprising considering the
uncertainties on the injection index are small. The CI model fitted to the
integrated spectrum is shown in Fig. 4.

Results of the fitting are summarised as follows (see also Tab. 2). Modelling
the integrated radio spectrum using the CI model gives a break frequency of
$\nu_{\mathrm{b}}=1.78_{-0.15}^{+0.17}$ GHz, and a remnant ratio of
$R_{\mathrm{rem}}=0.23\pm 0.02$. Using instead the TCI model, gives a
marginally different break frequency of
$\nu_{\mathrm{b}}=1.26_{-0.05}^{+0.06}$ GHz, and remarkably, a remnant ratio
of $R_{\mathrm{rem}}=0.23\pm 0.02$, consistent with the value constrained from
the standard CI model.

Table 2: The spectral modelling attributes of J2253-34, described in Sect. 3.4. Using the $8^{\prime\prime}$ radio images, the radio spectra arising from small regions across the lobes are jointly fit by the JP models, in order to optimise the injection index. The broadband integrated radio spectrum is then fit by the CI models, propagating the fitted injection index and uncertainties through a Gaussian prior, to estimate the break frequency and remnant ratio.  Model | Injection index | Break frequency | Remnant ratio
---|---|---|---
JP | $2.10\,\pm\,0.05$ | $-$ | $-$
TJP | $2.25\,\pm\,0.03$ | $-$ | $-$
CI | $-$ |  $1.78_{-0.15}^{+0.17}$ GHz | $0.23\pm 0.02$
TCI | $-$ |  $1.26_{-0.05}^{+0.06}$ GHz | $0.23\pm 0.02$

Finally, it is important to mention that the precise measurements of the break
frequency and remnant ratio of J2253-34 are only made possible by the
independently-fitted injection index. By inspecting the marginal distribution
of the fitted break frequency and remnant ratio, we find that their PDFs are
Gaussian-like, and importantly do not demonstrate any parameter degeneracy
worth consideration. These results are in stark contrast to those found by
blindly fitting each parameter of the CI models from the integrated radio lobe
spectrum. Here, we find that a single precise measurement of the injection
index can not be found; e.g. for the standard CI model, the fitted injection
index can plausibly range anywhere between $s\sim 2.1$ \- $2.3$. Consequently,
this introduces major uncertainties onto the fitted break frequency and
remnant ratio, which become largely degenerate with one another; e.g. for the
standard CI model, the break frequency demonstrates frequent local
minima/maxima between $\nu_{\mathrm{b}}\sim 2$ \- $5$ GHz, and the PDF of the
fitted remnant ratio is heavily asymmetric about its peak probable value of
$R_{\mathrm{rem}}\sim 0.38$. It is interesting that these results are found
even for an integrated radio spectrum as well-sampled as that of J2253-34, and
thus reinforce the approach described in this section.

##  4 DYNAMICAL MODELLING OF REMNANT LOBES

In the following section, we present a new dynamical model-based method that
measures the off-time in remnant lobes, and is capable of identifying
candidate remnant radio galaxies from observations. In Sect. 4.1, we present
the dynamical model that underpins this work. In Sect. 4.2, we justify why the
backflow of lobe plasma holds valuable constraints on the off-time, and
necessarily devise a tool to parameterise the surface brightness distribution
of radio lobes. In Sect. 4.3, we then outline the constraints necessary to
model the energetics of lobed radio-loud AGNs. Finally in Sect. 4.4, we
implement our new method to constrain the energetics of J2253-34, which we
then independently verify against the remnant ratio previously measured in
Sect. 3.

###  4.1 Radio AGN in semi-analytic environments

Throughout this work, the Radio AGN in Semi-analytic Environments (RAiSE;
Turner & Shabala, 2015; Turner et al., 2018a; Turner & Shabala, 2020) code is
used to model the dynamics of Fanaroff & Riley type-II (FR-II; Fanaroff &
Riley, 1974) radio lobes. In this work, we employ RAiSE both as a tool to
model the radio-loud AGN energetics, and to simulate mock radio source
populations.

####  4.1.1 The RAiSE model

The RAiSE code consolidates models describing the expansion of lobed FR-II/I
sources together with an observation-based treatment of the environments
within which the jets expand; RAiSE has been tested to show consistency with
hydrodynamical simulations, is able to reproduce the observed relationship
between radio luminosity, morphology and host galaxy properties, is consistent
with observed spectral ages maps of radio lobes, and is able to reconcile the
discrepancy between spectral and dynamical ages in powerful radio galaxies. In
this way, the temporal evolution of a radio source can be modelled for,
amongst other intrinsic properties, any combination of source age, jet kinetic
power, and intergalactic environment (forward modelling hereafter).

Quantifying an AGN environment requires knowledge of the distribution of gas
in the associated dark matter halo. Ideally this is done by constraining a gas
density and temperature profile through X-ray observations, however such
observations are largely unavailable across the majority of haloes, especially
at higher redshifts. Turner & Shabala (2015) showed that estimates of the halo
mass of the intergalactic environment can alternatively be inferred from the
optically-observed stellar content based on outputs of semi-analytic galaxy
evolution models (e.g. SAGE; Croton et al., 2006); however, the shape of the
gas density profile remains unconstrained by this method. RAiSE therefore
quantifies the intergalactic environment in one of two ways: 1) by directly
specifying a gas-density profile; or 2) by providing a mass of the dark matter
halo, distributing the total gas mass based on an observationally-based gas
density profile reported by Vikhlinin et al. (2006). In the latter example,
this can be done based on either the mean profile or by randomly sampling a
distribution of environment profiles from a galaxy formation model.

RAiSE implements a lobe emissivity model that considers radiative losses due
to synchrotron radiation and the inverse-Compton up-scattering of CMB photons,
as well as the adiabatic losses related to the expansion of the lobes. This
emissivity calculation is applied to small packets of electrons that are
shock-accelerated at different points in the evolutionary history of the lobe.
As a result, RAiSE considers the temporal evolution in the magnetic field;
i.e. from the time electrons are shock-accelerated to the time they release
their energy as synchrotron radiation. The magnetic field strength is derived
from the lobe dynamics by assuming a constant ratio between the energy
densities in the magnetic field and particles; referred to as the
equipartition factor.

Over the past decade, models describing these processes have been
incrementally added to the RAiSE code. Turner & Shabala (2015) provided the
original implementation of a radio source model within a semi-analytic
environment, and correspondingly modified the synchrotron emissivity model for
the lobes. Turner et al. (2018a) used hydrodynamical simulations to spatially
map the synchrotron emissivity across the lobes, allowing them to produce
surface brightness maps of simulated radio sources; their model implemented
tracer particles from hydrodynamical simulations to probe the average
locations and ages of injected particles. Turner (2018) extended the dynamical
and synchrotron emissivity models to account for the remnant phase during
which the jets are inactive. Of relevance to this work, the dynamics presented
by Turner & Shabala (2020) offer a slight modification to include the effect
of a shocked gas shell.

The latest version of the RAiSE code666https://github.com/rossjturner/RAiSEHD,
as used in this work, is informed by the dynamics of Turner & Shabala (2020),
and implements the method of Turner et al. (2018a) to produce intensity maps
based on tracer particles from the higher-resolution hydrodynamical simulation
of Yates-Jones et al. (2022). By using individual hydrodynamical simulation
particles rather than population averages, this enables RAiSE to synthesise
hydrodynamical-based, high-resolution maps of the radio lobes at any observing
frequency (e.g. see Fig. 5).

####  4.1.2 Constraining a RAiSE model

The RAiSE model is characterised by a variety of model parameters that
represent intrinsic parameters of a radio source. These parameters are
constrained in one of three ways: 1) through a direct measurement (e.g.
injection index); 2) a prior probability distribution based on measurements of
well-studied objects (e.g. ambient density profile); or 3) by way of a
parameter inversion based on simulated observable attributes (e.g. jet power
and age from size and luminosity).

RAiSE-based parameter inversions have been implemented by several previous
authors to constrain the energetics of radio-loud AGNs. Turner & Shabala
(2015) demonstrated the ability to recover the jet powers and lifetimes of
radio sources using only their observed physical size and monochromatic radio
luminosity; in their work, halo masses for each source were inferred from the
stellar masses of their hosts using SAGE, and intrinsic parameters such as the
equipartition factor and injection index were fixed to typical values. Turner
et al. (2018b) were additionally able to recover the lobe equipartition
factors by including the break frequency as an observational constraint.
Turner et al. (2020) demonstrated that using only a reference flux density,
angular size and width, break frequency and injection index, the cosmological
redshift of radio sources can be constrained from the present epoch to the
early-universe.

These results underline that the constraining power of the RAiSE-based
parameter inversion is primarily limited by the input set of observational
constraints. It is worth nothing that previous such applications of the RAiSE
model have used attributes measured from observations with poor resolution.
The capabilities of the current RAiSE model create a new opportunity to
constrain intrinsic radio source parameters using additional attributes
related to the observed surface brightness distribution of radio lobes. In the
following sections, we exploit this opportunity to explore a new method for
modelling the energetics of remnant AGNs.

Figure 5: Intensity maps predicted by RAiSE for a redshift $z=0.1$ radio
source with a jet-power of $Q=10^{38}$ W, source age of $\tau=100$ Myr, and a
remnant ratio of $R_{\mathrm{rem}}=0.25$. The radio source is shown for 150
MHz (left panels) and 1.4 GHz (right panels), each of which is shown for a
native resolution corresponding to $1.024\times 10^{8}$ particles (upper
panels) and for a convolved 10′′ image (lower panels). The pixel scale in the
convolved image is $1^{\prime\prime}$. Following the method outlined in Sect.
4.2.2, the observed surface brightness distribution is optimally fit by a
skewed two-dimensional Gaussian; this is shown by a solid black contour, which
represents the edge bounding 95% of the Gaussian volume. Images of each fitted
Gaussian are truncated at the core, resulting in the apparent abrupt cut-off
at 150 MHz.

###  4.2 The surface brightness distribution of radio lobes

In the absence of detailed spectral modelling, alternative constraints are
needed to aid the modelling of the energetics of remnant AGNs. Here, we
propose that the backflow of plasma in radio galaxy lobes offers such
constraints, and subsequently develop a tool to parameterise the observed
surface brightness distribution of radio lobes.

####  4.2.1 Constraints embedded in the lobe backflow

The plasma contained within FR-II radio lobes is injected throughout the
active phase at the hotspot, and is carried away from the hotspot by the
backflow (e.g. Turner et al., 2018a). This leads to a gradient in age across
the lobes, generally in the form of an approximately linear increase in
spectral age away from the hotspot. The observed shape of the emitting regions
will vary depending on observing frequency; the low-frequency radiation
produced by older electrons results in wider emitting regions along the jet-
axis and transverse-axis, whereas at higher frequencies the bias towards
observing younger electrons results in narrower emitting regions localised
near the injection site. This is demonstrated in Fig. 5. Importantly, the
manner by which the observing frequency modifies the observed spatial
distribution of synchrotron-emitting electrons, should depend on factors such
as the plasma age distribution and lobe magnetic field strength. In the
absence of spectral modelling, it is therefore reasonable to suggest that
well-resolved observations of remnant radio lobes should offer the constraints
necessary to help disentangle their active and remnant timescales.

####  4.2.2 Skewed two-dimensional Gaussian lobes

We require attributes to describe the observed surface brightness distribution
of radio lobes, which can be compared between observations and simulations. We
note that the reproduction of fine details is neither important nor captured
by the RAiSE simulations (e.g. knots in the jets, or slight bending of the
jets). Implementing a skewed Gaussian in two spatial dimensions is appropriate
for the following reasons: 1) a Gaussian is able to describe the fall-off in
brightness as a function of position; 2) in two dimensions, a Gaussian is able
to characterise a measure of length of the lobe along the jet-axis, as well as
a width along the transverse-axis; 3) introducing a ‘skew’ term along the
major axis handles the asymmetry in brightness along the jet-axis (e.g.
emission is typically brightest near the hotspot and falls off asymmetrically
in each direction). In this way, the vastly different morphologies of, for
example, a high-powered radio source (with a sharp fall-off in brightness away
from the hotspot) and a low-powered aged remnant (with relatively smooth
surface brightness profiles) can be characterised in a manner that is
consistent.

To fit the skewed Gaussian to the observed surface brightness distribution of
a radio lobe, the following free parameters are necessary:

  * •

$S_{0}$, the amplitude of the Gaussian;

  * •

$\mu_{x}$ & $\mu_{y}$, the coordinates of the peak amplitude in each
dimension;

  * •

$\sigma_{x}$ & $\sigma_{y}$, the standard deviation in each dimension;

  * •

$\theta$, the major axis position angle;

  * •

$\beta_{x}$, the skew term along the Gaussian major axis.

Since RAiSE simulates a jet in an known orientation and positions the SMBH at
the center of the image, the angle of rotation ($\theta=0$) and $\mu_{y}$ are
set by default (and therefore do not need solving). A simple least-squares
algorithm is used to optimise the fit by minimising the residuals between the
radio lobe surface brightness distribution and Gaussian model. For stability,
regression is performed to the square-root of the input image; this ensures
fitting to the general shape, and down-weights bias towards the brightest
pixels. The resulting Gaussian fit can then be used to parameterise key
attributes related to the observed surface brightness distribution of a radio
lobe (see Sect. 4.3.1 below). We make this code publicly available on
GitHub777https://github.com/benjaminquici/skewed-gaussian-lobes.

###  4.3 Radio source attributes

In this work, radio source attributes refer to the observable quantities that
describe either the spatial, photometric, or spectral properties of radio
galaxy lobes. By implementing these attributes in a RAiSE-based parameter
inversion, they facilitate the comparison between observed radio sources and
those simulated by the RAiSE model. However, images of real radio sources are
‘degraded’ by the sensitivity and resolution of the observations; degrading
simulated images with these same observing limitations therefore enables a
more meaningful comparison. In the following section, we describe the methods
for measuring these attributes from observed and simulated images.

Figure 6: Schematic showing the key quantities necessary for computing the
spatial attributes described in Sect. 4.3.1. The semi-transparent purple
emission represents the right lobe of the 1.4 GHz convolved image from Fig 5.
A skewed two-dimensional Gaussian is fit to the lobe surface brightness
distribution following the method outlined in Sect. 4.2.2, where $S_{0}$ gives
the amplitude of the fit. A solid black contour denotes the fitted Gaussian at
a critical brightness level of $S_{0}e^{-2}$, and bounds 95% of the Gaussian
volume.

####  4.3.1 Spatial attributes

Spatial attributes are measured based on the Gaussian fitting method described
in Sect. 4.2.2. We do this by first considering the boundary of the 2D ellipse
formed by slicing the fitted Gaussian at some critical brightness level (e.g.
see Fig 5). We define this critical level as $S_{0}e^{-2}$, above which 95% of
the Gaussian volume is contained. We stress that the arbitrary selection of
this value does not matter; instead, what is important is that the same value
is used when comparing both observed and simulated radio images, such that the
resulting attributes are measured consistently. The following distances are
measured by considering the boundary of the 2D ellipse (Fig. 6): the length of
the lobe, $D_{l}$; the visible length of the lobe, $\bar{D_{l}}$; and the
radius of the lobe, $R_{l}$. We consider the backflow in each lobe separately,
and thus bound the visible lobe length as $0<\bar{D_{l}}\leq D_{l}$. From
here, we define four key spatial attributes, including the linear size, $D$,
and observed axis ratio, $A$, as shown in Tab. 3.

Importantly, to constrain the duration of the remnant phase we develop two new
measures of “extent”, which parameterise the observed backflow of lobe plasma.
The “extent 1” attribute, $E$, considers what fraction of the total lobe
length, parallel to the jet axis, is visible; e.g. how far back the lobe
visibly sweeps towards the core. Lobes with emission localised near the
hotspot will measure $E\rightarrow 0$, which is expected if a majority of the
lobe plasma energy is stored in the most freshly injected electrons.
Conversely, $E\rightarrow 1$ implies the observed radio emission sweeps all
the way to (and beyond) the core, which is expected if a majority of the lobe
plasma energy is stored in the older synchrotron-emitting electrons. The
“extent 2” attribute, $\bar{E}$, considers how bright the lobes are at the
core, relative to their peak brightness. We expect this attribute to be
particularly useful at disentangling sources with backflow observed all the
way to the core, where the measurement of the “extent 1” attribute will likely
just reflect its upper limit of $E=1$.

Table 3: Spatial attributes evaluated using the quantities shown in Fig. 6 Spatial attribute | Symbol | Expression | Limit
---|---|---|---
linear size | $D$ | $D_{l}^{\mathrm{{lobe1}}}+D_{l}^{\mathrm{{lobe2}}}$ | $-$
observed axis ratio | $A$ | $D_{l}/R_{l}$ | $\geq 1$
extent 1 | $E$ | $D_{l}/\bar{D_{l}}$ |  $0$ \- $1$
extent 2 | $\bar{E}$ | $S_{\mathrm{core}}/S_{0}$ |  $0$ \- $1$

####  4.3.2 Photometric attributes

A total radio luminosity, $L_{\nu}$, offers a powerful constraint on radio
source energetics. For an observed radio source, the total radio luminosity
can be derived from the total radio flux density by integrating over a region
of the source bounded by a minimum brightness criterion; e.g. the region of
the source above some signal to noise threshold. This means that, for a fixed
sensitivity, the measured integrated flux density will depend on the surface
brightness distribution of the lobes. The surface brightness maps produced by
the RAiSE model allow us to include the effect of instrumental noise in the
computation of the integrated flux density, thus accounting for this
systematic bias. This is done by populating the simulated surface brightness
maps with a Gaussian noise term, then the total observed radio luminosity is
calculated considering only pixels brighter than $3\sigma$ (where $\sigma$
represents the level of Gaussian noise). This process of adding the noise is
repeated over 16 iterations to sample a mean integrated flux density, and to
characterise an uncertainty, $dL_{\nu}$, based on the standard deviation
across each evaluation.

By measuring the integrated flux density over two observing frequencies,
$\nu_{0}$ and $\nu_{1}$, a two-point spectral index can be calculated through:

$\alpha~{}=~{}\log(L_{\nu,1}/L_{\nu,0})/\log(\nu_{1}/\nu_{0}),$ (10)

which has an uncertainty corresponding to:

$d\alpha~{}=~{}\big{(}1/\ln({\nu_{0}/\nu_{1}})\big{)}~{}\sqrt{(dL_{\nu,0}/L_{\nu,0})^{2}+(dL_{\nu,1}/L_{\nu,1})^{2}}.$
(11)

With at least three observing frequencies a spectral curvature, SPC, can be
calculated through: $\mathrm{SPC}=\alpha_{0}-\alpha_{1}$, which has an
uncertainty corresponding to
$d\mathrm{SPC}~{}=~{}\sqrt{(d\alpha_{0})^{2}+(d\alpha_{1})^{2}}$.

Table 4: The measured radio source attributes of J2253-34, which are measured following the methods outlined in Sect. 4.3. In Sect. 4.4.2 these attributes are compared with corresponding RAiSE model outputs, using the scale shown in column 4, to constrain the energetics via a parameter inversion. To explore the constraints embedded in the lobe backflow, attributes are compared in four different combinations (columns 5-8). Tick ($\checkmark$) and cross (✗) markers represent whether the attributes are included or excluded from the fitting, respectively. Radio source attribute | Symbol | Measurement | Scale | full fit | spectral | spatial | in band
---|---|---|---|---|---|---|---
radio luminosity | $L_{1.4}$ |  $(6.48\pm 0.25)\times 10^{24}$ W Hz-1 | $\log_{10}$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$
${}^{\text{(measured at 1.42\,\text{GHz})}}$ |  |  |  |  |  |  |
largest linear size | $D_{1.4}$ | 326 kpc | $\log_{10}$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$
${}^{\text{(measured at 1.42\,\text{GHz})}}$ |  |  |  |  |  |  |
break frequency | $\nu_{\mathrm{b}}$ |  $1.26^{+0.06}_{-0.05}$ GHz | $\log_{10}$ | $\checkmark$ | $\checkmark$ | ✗ | ✗
observed axis ratio |  $A_{1.4}$, $A_{5.5}$ |  $2.4\pm 0.1$, $2.7\pm 0.1$ | linear | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$
${}^{(1.42,~{}5.50\,\text{GHz})}$ |  |  |  |  |  |  |
extent 1 |  $E_{1.4}$, $E_{5.5}$ |  $0.9\pm 0.05$, $0.8\pm 0.05$ | linear | $\checkmark$ | ✗ | $\checkmark$ | ✗
${}^{(1.42,~{}5.50\,\text{GHz})}$ |  |  |  |  |  |  |
extent 2 |  $\bar{E}_{1.4}$, $\bar{E}_{5.5}$ |  $0.76\pm 0.05$, $0.65\pm 0.05$ | linear | $\checkmark$ | ✗ | $\checkmark$ | ✗
${}^{(1.42,~{}5.50\,\text{GHz})}$ |  |  |  |  |  |  |
spectral index | $\alpha_{150}^{1420}$ | $0.64\pm 0.03$ | linear | ✗ | ✗ | ✗ | $\checkmark$
${}^{(150-1420)\,\text{MHz}}$ |  |  |  |  |  |  |
spectral curvature | $\mathrm{SPC}_{150}^{5500}$ | $0.53\pm 0.04$ | linear | ✗ | ✗ | ✗ | $\checkmark$
${}^{(150-1420-5500)\,\text{MHz}}$ |  |  |  |  |  |  |

####  4.3.3 Spectral modelling attributes

The injection index and remnant ratio attributes, derived from observations by
modelling the spectra of remnant radio galaxy lobes (Sect. 3), represent
intrinsic parameters of the RAiSE model. On the other hand, a break frequency
can be predicted by the RAiSE model, which offers valuable constraints on the
lobe equipartition factor. While it is possible to measure a RAiSE-predicted
break frequency by modelling a simulated radio spectrum, the computational
cost of simulating a radio source at many different frequencies makes this
method inefficient. For any given set of model parameters, the RAiSE code uses
the lobe pressure and equipartition factor to compute the lobe magnetic field
strength. In turn, the field strength is used through Eqn. 3 to calculate the
corresponding RAiSE-predicted break frequency, where the redshift and source
age are those used as RAiSE model inputs. By modelling the spectra of RAiSE-
simulated radio sources, Turner et al. (2018b) showed that the break frequency
computed from the magnetic field strength is consistent with that derived by
modelling their synthetic radio spectra (e.g. see their Fig. 3).

###  4.4 Modelling the energetics of J2253-34

We are now equipped with the tools needed to model the energetics of J2253-34.
To begin with, we model the energetics closely following the dynamical model
based method presented by Turner (2018). Here, the remnant ratio is already
known (e.g. from Sect. 3.4.2), whilst the jet power, source age, and
equipartition factor are principally constrained from a radio luminosity,
physical size, and break frequency. We compare this established technique to
our proposed method using the spatial attributes of the lobes (namely extent)
to additionally constrain the remnant ratio. Finally, we provide our most
precise constraints on the energetics by combining both the spectral and
spatial attributes to best understand the physics of J2253-34.

####  4.4.1 Observed properties of J2253-34

The radio galaxy J2253-34 was classified as a remnant by Quici et al. (2021)
based on the absence of radio emission from the core, and further supported by
the observed steepening in the integrated radio spectrum. In that work, the
radio source was unambiguously associated with a $z=0.2133$ host galaxy, for
which a halo mass of $M_{\mathrm{H}}=10^{13.5}$ M⊙ was inferred based on the
stellar mass of the host.

The radio source attributes of J2253-34 are measured following the methods
outlined in Sect. 4.3, and are summarised in Tab. 4. All spatial attributes
are measured from the $8^{\prime\prime}$ images presented in Sect. 2, and we
use the injection index, break frequency, and remnant ratio fitted in Sect.
3.4 by the Tribble spectral ageing models (e.g. row 2 & 4 of Tab. 2). In the
following section, we implement these attributes to constrain the energetics
of J2253-34 via a parameter inversion.

Table 5: The intrinsic parameters of J2253-34, constrained via the four
separate methods outlined in Tab. 4.
† parameterised by the CI model. Param. | full fit | spectral | spatial | in band
---|---|---|---|---
$Q$ | $3.98^{+0.36}_{-0.40}$ | $3.98^{+0.36}_{-0.40}$ | $6.32^{+2.11}_{-1.63}$ | $1.0$ \- $7.9$
${}^{(\times 10^{38}\text{ W})}$
$\tau$ | $56.2^{+2.62}_{-1.28}$ | $56.2^{+2.62}_{-1.28}$ | $50.1^{+4.2}_{-3.8}$ | $31$ \- $158$
${}^{(\text{Myr})}$
$B$ | $1.98^{+0.16}_{-0.15}$ | $1.98^{+0.16}_{-0.15}$ | $1.41^{+0.30}_{-0.25}$ | $2$ \- $14$
${}^{(\times 10^{-10}\text{ T})}$
$R_{\mathrm{rem}}$ | ${}^{\dagger}0.23^{+0.02}_{-0.02}$ | ${}^{\dagger}0.23^{+0.02}_{-0.02}$ | $0.26^{+0.04}_{-0.03}$ |  $0.1$ \- $0.3$

####  4.4.2 RAiSE-based parameter inversion of J2253-34

RAiSE can model the temporal evolution of radio lobes for any set of intrinsic
radio source parameters, and subsequently generate mock images to make
predictions for their respective radio source attributes. The intrinsic
parameters of any radio source can therefore be estimated via a RAiSE-based
parameter inversion; observed radio source attributes are compared to a
corresponding set of mock attributes generated over a multi-dimensional set of
intrinsic parameters (e.g as is done by; Turner & Shabala, 2015; Turner et
al., 2018b; Turner, 2018; Shabala et al., 2020).

For a fixed redshift of $z=0.2133$ and halo mass of $M_{\mathrm{H}}=10^{13.5}$
M⊙, a grid of RAiSE models are created from the two-sided jet kinetic power,
$Q$; the dynamical source age, $\tau$; the equipartition factor, $q$ (which
parameterises the ratio of energy densities in the magnetic field and
particles; Turner et al., 2018a); the initial supersonic phase value of the
axis ratio, $A_{i}$; the remnant ratio, $R_{\mathrm{rem}}$; and the injection
index, $s$. We use the results of Turner et al. (2018b), who fit the intrinsic
properties of 3C radio sources, to ensure our model grid covers plausible
parameter values. In this work, RAiSE models are evolved between 10 Myr to 1
Gyr, recording outputs at every $\Delta(\log\tau)=0.05$ dex. We uniformly
sample $\log Q$ between $[37,40]$ with a precision of $0.05$ dex; $\log q$
between $[-2.7,-0.7]$ with a precision of $0.05$ dex; $A_{i}$ between $[2,8]$
with a precision of $0.25$; $R_{\mathrm{rem}}$ between $[0,5]$ with a
precision of $0.01$; and $s$ between $[2.01,2.61]$ with a precision of $0.1$.

To explore the constraining power of the extent attributes, we perform the
parameter inversion based on four separate combinations of the radio source
attributes, summarised in cols. $5$ \- $8$ of Tab. 4. The spatial method tests
a scenario where detailed spectral modelling (e.g. Sect. 3.4) can not be
performed; here, each intrinsic radio source parameter is explicitly fit via
the RAiSE-based parameter inversion, and the break frequency is omitted as a
constraint. To independently verify results of this method, the fitted
energetics are compared with that fitted by the full fit and spectral methods;
here, the injection index and, importantly the remnant ratio, are fixed to
known values, and the strong constraint from the break frequency is included.
The in band method considers a limiting case where no constraints on the
duration of the remnant phase are available.

To conduct the parameter inversion (for any of the above methods), simulated
images are used to measure the same attributes as those measured for J2253-34,
ensuring that the images are degraded by the same set of observing
constraints. Attributes are compared with model outputs through Eqn. 9, and we
follow the same process described in Sect. 3.2 to select the best-fit model,
and to estimate the uncertainty in each parameter.

The values fitted for the intrinsic parameters using each of these methods are
presented in Tab. 5. Consistent with the results of Turner (2018), we find
that the jet-power/age/equipartition-factor triplet, fitted by the spectral
method, is tightly constrained by the observed radio-luminosity/linear-
size/break-frequency. We inspect the marginal distributions of each fitted
parameter to find approximately Gaussian-like PDFs centered about the peak
probable value; i.e. the fitting appears stable and not degenerate. By
including constraints from the extents, we find that the energetics
constrained by the full fit method do not change. This suggests that the break
frequency is the dominant observable in this fitting process.

Using the source age constrained by full fit method, together with the remnant
ratio fit previously by the TCI model, the duration of the active and remnant
phases in J2253-34 are calculated as $t_{\mathrm{on}}=43.3\pm 4.3$ Myr and
$t_{\mathrm{rem}}=12.9\pm 1.2$ Myr, respectively. Importantly, through this
method we demonstrate the ability to accurately measure the duration of the
active phase, and to place a weak upper bound on the AGN jet duty cycle
($\delta\lesssim 1-R_{\mathrm{rem}}=0.77$).

Remarkably, we find that the energetics modelled by the spatial method are
consistent with those derived by the previous methods. Uncertainties in each
parameter are comparatively greater, as expected considering the constraints
are relaxed and an additional free parameter is solved for, however inspecting
the marginal distributions reveals similar Gaussian-like PDFs. What is
particularly of interest is that the dynamically-modelled remnant ratio,
$R_{\mathrm{rem}}=0.26\pm 0.02$, is consistent with that parameterised by the
TCI model within a level of $1\sigma$. Together with a source age of
$\tau=50.1^{+4.2}_{-3.8}$ Myr, the spatial method therefore estimates the
duration of the active and remnant phases as $t_{\mathrm{on}}=37.2\pm 4.3$ Myr
and $t_{\mathrm{rem}}=13.1\pm 1.5$ Myr, respectively. These timescales are
consistent with those fitted by the full fit method, suggesting that the
surface brightness distribution of the lobe backflow can in fact be used to
constrain the dynamics of the remnant phase.

Finally, we find that the stability of the model fitting breaks down severely
with the in band method. Inspecting the marginal distributions reveals a large
degeneracy in the fitted parameter space, where multiple unique values of
equal likelihood are found over a wide range of parameter values. We quote
these ranges in Tab. 5, by considering their minimum/maximum values. This
tells us that the spectral index and curvature, alone, do not offer the
constraining power needed to tightly constrain the energetics of J2253-34.

##  5 MOCK RADIO SOURCE POPULATIONS

Figure 7: Results of the parameter inversion of the $z=0.1$ (upper row) and
$z=0.5$ (lower row) mock catalogues are shown. Each plot demonstrates the
fitted (vertical axis) versus mock (horisontal axis) parameter for the jet
power (first column), source age (second column), equipartition factor (third
column), and remnant ratio (fourth column). Each marker represents a unique
mock radio source, where unfilled markers represent sources fainter than
$S_{150}=100$ mJy or smaller than $\theta=60^{\prime\prime}$. A solid black
line is shown to represent a 1:1 relation between the fitted and mock
quantities.

The ability to measure the off-time in remnant radio galaxies, using
attributes measured exclusively below $\sim 1.4$ GHz, is particularly
attractive considering the parameter space unlocked by wide-area sky surveys
such as the LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al., 2019, 2022),
the MeerKAT International GHz Tiered Extragalactic Exploration (MIGHTEE;
Jarvis et al., 2016) survey, the APERture Tile In Focus (Apertif Oosterloo et
al., 2010) survey, and EMU. However, the generalisability of this new method
is currently unclear, as we do not know whether the constraining power of the
spatial attributes, measured using telescope-degraded images, deteriorates for
certain intrinsic source parameters; e.g. lower jet powers, older sources, or
denser environments. In Sect. 5.1, we use our method to perform a parameter
inversion of mock radio source populations simulated at 150 MHz and 1.4 GHz.
In Sect. 5.2, we compare the input and recovered parameters to examine the
internal consistency of our method within the RAiSE model; e.g. does the
method recover the input parameters. In Sect. 5.3 we then examine whether the
remnant ratio fitted by this method can be used to confidently select between
active and remnant radio lobes.

###  5.1 Constructing mock radio source catalogues

For this analysis, mock radio source populations are simulated using the RAiSE
model. We allow the following set of intrinsic parameters to vary over a wide
range of values, typically associated with FR-IIs:

  * •

_jet power:_ $Q\in[10^{37},~{}10^{40}]\;\rm W$, to represent those typically
associated with FR-II jets (Turner & Shabala, 2015)

  * •

_source age:_ $\tau\in[10^{7},~{}10^{8.5}]\;\rm yr$. Younger sources will
likely fall below extended source criteria (see Fig. 7), and a genuine dearth
of old ($>400$ Myr) radio sources has been reported by Shabala et al. (2020)
in LOFAR data.

  * •

_remnant ratio:_ $R_{\mathrm{rem}}\in[0,~{}0.5]$, consistent with the
observational bias towards selecting recently switched-off remnants (Mahatma
et al., 2018; Turner, 2018; Jurlin et al., 2020; Shabala et al., 2020; Quici
et al., 2021).

  * •

_equipartition factor:_ $q\in[10^{-2.7},~{}10^{-0.7}]$, which is approximately
Gaussian-distributed about $\log q\approx-1.7$ (e.g. see; Turner et al.,
2018b).

  * •

_initial axis ratio:_ A${}_{i}\in[2.5,~{}4.5]$, reflecting a subset of typical
values for the 3C sample (e.g. see; Turner & Shabala, 2015).

  * •

_injection index:_ $s\in[2.01,~{}2.61]$, where the lower limit comes from
predictions of first-order Fermi acceleration (Kardashev, 1962; Pacholczyk,
1970), and the upper limit from $\alpha_{\mathrm{inj}}$ not being too steep.

Figure 8: The fitted (vertical axes) versus input (horisontal axes) remnant
ratio of the $z=0.1$ (upper panel) and $z=0.5$ (lower panel) mock catalogues.
Each column represents a log-uniform source age bin. Gray regions represent a
$1\sigma$ uncertainty in the scatter, corresponding to the unflagged sources
(filled markers). A tightening in the scatter is seen as the source age
increases, likely reflecting the greater magnitude of radiative losses
discussed in Sect. 5.2. The dearth of sources with higher values of
$R_{\mathrm{rem}}$ in the $z=0.5$ sample arises due to the complete depletion
of electrons capable of producing emission at these frequencies (e.g. see Fig.
4 of Turner 2018).

Two separate mock catalogues are created at redshift $z=0.1,0.5$ (hereafter
the $z=0.1$ and $z=0.5$ catalogues). For each catalogue, 5000 mock radio
sources are forward modelled using pseudo-random values for the intrinsic
parameters. We use a representative $M_{\mathrm{H}}=10^{13}$ M⊙ halo mass
(e.g. a typical environment of an FR-II, Turner & Shabala 2015), and simulate
a random gas-density profile based on a prior distribution of environment
profiles (as per discussion in Sect. 4.1.1).

A corresponding grid of RAiSE models is simulated with known values for each
intrinsic parameter, so that a parameter inversion can be performed on each
mock radio source catalogue. Our model grids are bound by the parameter space
limits of the mock sample, however with uniform spacings as follows:
$\Delta(\log Q)=0.05$ dex; $\Delta(\log\tau)=0.025$ dex; $\Delta
R_{\mathrm{rem}}=0.01$; $\Delta(\log q)=0.05$ dex; $\Delta A_{i}=0.25$; and
$\Delta s=0.1$. Sources are simulated for the same halo mass, however here we
take the mean of the environment profiles.

Bayesian parameter estimation is conducted as follows. Images are synthesised
at 150 MHz and 1.4 GHz, and degraded by a set of survey limitations shown in
Tab. 6. This is done both for mock radio sources and those stored in the grid
as model outputs. Telescope-degraded images are then used to measure the radio
source attributes shown in Tab. 7. Following the parameter inversion described
in Sect. 4.4.2, the attributes measured for mock radio sources are compared
with those used as model outputs, in order to constrain their intrinsic
parameters. Here, we simultaneously fit for the pseudo-randomly sampled
parameters, and treat the redshift and halo mass as known.

Table 6: The set of observing conditions used to sample mock radio source populations simulated in Sect. 5.1.  Survey type | Freq. | Sensitivity | Resolution
---|---|---|---
(MHz) | ($\mu$Jy beam-1) | (′′)
LOFAR-like | 150 | 70 | 6
${}^{\emph{(wide-band)}}$ |  |  |
LOFAR-like | 130 | 100 | 6
${}^{\emph{(sub-band)}}$ | 170
MeerKAT-like | 1400 | 6 | 10

Sources are flagged if their angular size falls below
$\theta=60^{\prime\prime}$, or if their 150 MHz integrated flux densities fall
below $S_{150}=100$ mJy. This decision is motivated by the typical current
observational constraints used to construct complete radio galaxy samples
(e.g; Brienza et al., 2017; Godfrey et al., 2017; Mahatma et al., 2018; Jurlin
et al., 2020; Quici et al., 2021). The results of this parameter inversion are
explored in the following sections.

Finally, we note that the Karl G. Jansky Very Large Array Sky Survey (VLASS;
Lacy et al., 2020) represents a high-frequency ($2-4$ GHz), high-resolution
($2.^{\prime\prime}5$), all-sky (decl.$>-40^{\circ}$) radio survey, which is
expected to achieve $\sigma=70\mu$Jy beam-1 sensitivity. Despite marginally
increasing the frequency leverage over 1.4 GHz, we do not consider a VLASS-
like observation due to the comparatively poorer surface brightness
sensitivity to extended radio emission; without a handle over the fainter
(older) lobe electrons, we expect the observed backflow to be much less
constraining on the remnant off-time.

###  5.2 Fitted energetics

We expect that the parameter inversion may be less effective for some sets of
parameters due to the non-linear application of radiative losses across our
parameter space. In particular, the two extent attributes, and to a lesser
degree the lobe width attribute, rely on the surface brightness distribution
of ageing synchrotron-emitting lobe plasma. The ageing is quantified through
the time-derivative of the Lorentz factor of a given packet of synchrotron-
emitting electrons, $\gamma$, as follows (Kaiser et al., 1997):

$\frac{d\gamma}{dt}=-\frac{a_{1}\gamma}{3t}-\frac{4\sigma_{\mathrm{T}}\gamma^{2}\,(u_{B}+u_{C})}{3m_{e}c},$
(12)

where the first and second terms encode adiabatic and radiative losses
respectively. Here, the derivative is taken with respect to the age of the
injected plasma, $t$, where $a_{1}$ gives the exponent with which the volume
of the lobe grows adiabatically, $\sigma_{\mathrm{T}}$ is the cross section
for Thompson scattering, and $u_{B}$ & $u_{C}$ are the energy densities of the
magnetic field and CMB respectively.

Table 7: The set of radio source attributes used to constrain the energetics of the mock radio source populations described in Sect. 5.1. Attribute | Scale
---|---
150 MHz radio luminosity | $\log_{10}$
largest linear size | $\log_{10}$
observed axis ratio (150, 1400) MHz | linear
extent 1 (150, 1400) MHz | linear
extent 2 (150, 1400) MHz | linear
($150-1420$) MHz spectral index | linear
($150-1420$) MHz spectral curvature | linear

Comparisons between their mock and fitted intrinsic parameters are presented
in Fig. 7 for the two mock catalogues described in Sect. 5.1. Firstly,
although with varying degrees of confidence (as indicated by the scatter), we
find that the fit of each parameter broadly closely matches the 1:1 relation
with its mock equivalent. The jet power and source age are fit with greater
confidence than the equipartition factor and remnant ratio, demonstrating the
former two property’s dependence on observable parameters; e.g. the jet power
and source age are tightly constrained by the radio luminosity and linear size
respectively, at least for active lobed sources (Turner & Shabala, 2015).
However, it is encouraging to find that the equipartition factor and remnant
ratio are generally recovered as well.

The jet power is less confidently predicted using the RAiSE-based parameter
inversion for lower-powered sources, as shown in Fig. 7. This is unlikely a
selection effect, in which the stochastic perturbations due to instrumental
noise are more significant for fainter sources associated with lower-powered
jets, since the same sizeable broadening is observed in the $z=0.1$ (brighter)
catalogue. Instead, this is likely due to the strong correlation between jet
power and the lobe magnetic field strength, and thus energy density $u_{B}$
(e.g. see Eqn. A5 of Turner et al. 2018a). The rate of energy loss from the
synchrotron-emitting electrons is therefore expected to be much greater in
higher jet power sources (Eqn. 12), thus leading to a greater rate of change
in the extent attributes as the remnant ages. Qualitatively this is consistent
with the results of Turner (2018), who find a rapid dimming in higher-powered
remnants in stark contrast to lower-powered sources whose flux density is
maintained over a greater duration (see their Fig. 4).

We also find that the energetics of the $z=0.5$ sample are fit with greater
confidence (see also discussion in next paragraph). Again referring to Eqn.
12, this result is unsurprising; the energy density of the CMB scales with
redshift as $u_{c}\mathchoice{\mathrel{\raise 1.07639pt\hbox{\hbox
to0.0pt{\hbox{$\displaystyle\propto$}\hss}\lower
4.03563pt\hbox{$\displaystyle\sim$}}}}{\mathrel{\raise 1.07639pt\hbox{\hbox
to0.0pt{\hbox{$\textstyle\propto$}\hss}\lower
4.03563pt\hbox{$\textstyle\sim$}}}}{\mathrel{\raise 0.75346pt\hbox{\hbox
to0.0pt{\hbox{$\scriptstyle\propto$}\hss}\lower
2.82494pt\hbox{$\scriptstyle\sim$}}}}{\mathrel{\raise 0.5382pt\hbox{\hbox
to0.0pt{\hbox{$\scriptscriptstyle\propto$}\hss}\lower
2.0178pt\hbox{$\scriptscriptstyle\sim$}}}}z^{4}$, meaning the magnitude of the
energy losses sustained throughout the remnant phase will be greater.

We investigate if the sizeable scatter in the remnant ratio is partially due
to one of the other model parameters. In particular, we expect the radiative
losses in the remnant phase will be related to the off-time, rather fractional
time spent in the remnant phase (e.g. integrating Eqn. 12 for constant volume;
i.e. $a_{1}=0$). As a result, we expect the synchrotron losses to reduce the
extent attributes in approximate proportion to the time spent in the off-
phase. The ability of the RAiSE-based parameter inversion to constrain the
remnant ratio is therefore expected to correspondingly be much less effective
in young sources (i.e. the off-time is smaller for a given remnant fraction).
To consider this potential sensitivity in our method, we group the mock
catalogues into four source age bins uniformly-spaced in $\log\tau$. The
fitted remnant ratio, in each source age bin, compared to the mock remnant
ratio is presented in Fig. 8. We find that the scatter in the predicted
remnant ratio reduces considerably for the older source age bins, as expected
from the above discussion.

Figure 9: Normalised count representing the distribution of fitted remnant
ratios of mock samples with $R_{\mathrm{rem}}=0$ (black),
$R_{\mathrm{rem}}=0.1$ (purple), and $R_{\mathrm{rem}}=0.2$ (blue). We set the
bin size in each distribution to $\Delta R_{\mathrm{rem}}=0.01$, to reflect
spacings in the grid. Each sample is simulated at redshift $z=0.1$, and only
the oldest source-age bin ($133~{}\leq~{}\tau~{}\leq~{}316$) Myr is shown.
Solid and dotted vertical lines are shown to represent the 68th and 95th
percentiles in each sample, respectively; these represent upper bounds for the
active population, and lower bounds for the remnant samples.

Overall, the results presented here show promise for modelling the energetics
of remnant AGNs, and notably, using spatial attributes of their lobes to
constrain the remnant phase. The predictive power of this technique improves
when the magnitude of the energy losses sustained during the remnant phase are
greater; this is unsurprising, as the technique relies on the reduction in
extent to constrain the ageing during this phase. The method makes reliable
estimates of the remnant ratio associated with older sources, in contrast to
young sources for which the fitted remnant ratio is less trustworthy.

###  5.3 Identifying remnant radio galaxies

Table 8: The recovery rate of active and remnant radio sources simulated in Sect. 5.3. Three mock radio source catalogues ($N=5000$) are created at $z=0.1$ and $z=0.5$, with remnant ratios fixed to $R_{\mathrm{rem}}=0,0.1,0.2$, respectively. For each redshift (column 1) and source age bin (column 2), the $68^{\mathrm{th}}$ and $95^{\mathrm{th}}$ percentiles of the fitted remnant ratio distribution are reported (columns $3-8$). Mock sample |  | Active | Remnant $(R_{\mathrm{rem}}=0.1)$ | Remnant $(R_{\mathrm{rem}}=0.2)$
---|---|---|---|---
Redshift | Source age (Myr) |  | (68th perc.) | (95th perc.) | (68th perc.) | (95th perc.) | (68th perc.) | (95th perc.)
$z=0.1$ | $10.0\leq\tau<23.7$ |  | $\leq 0.00$ | $<0.05$ | $>0.07$ | $>0.02$ | $>0.15$ | $>0.08$
$23.7\leq\tau<56.2$ |  | $\leq 0.00$ | $<0.12$ | $>0.07$ | $>0.01$ | $>0.15$ | $>0.06$
$56.2\leq\tau<133.0$ |  | $<0.01$ | $<0.10$ | $>0.07$ | $>0.01$ | $>0.16$ | $>0.06$
$133.0\leq\tau\leq 316.0$ |  | $<0.01$ | $<0.05$ | $>0.07$ | $>0.03$ | $>0.18$ | $>0.12$
$z=0.5$ | $10.0\leq\tau<23.7$ |  | $\leq 0.00$ | $<0.05$ | $>0.06$ | $>0.01$ | $>0.14$ | $>0.03$
$23.7\leq\tau<56.0$ |  | $\leq 0.00$ | $<0.03$ | $>0.08$ | $>0.05$ | $>0.17$ | $>0.13$
$56.2\leq\tau<133.0$ |  | $<0.01$ | $<0.03$ | $>0.09$ | $>0.06$ | $>0.19$ | $>0.17$
$133.0\leq\tau\leq 316.0$ |  | $\leq 0.00$ | $<0.02$ | $>0.09$ | $>0.08$ | $>0.19$ | $>0.18$

As we have demonstrated the ability of our method to estimate the remnant
ratio with reasonable accuracy, this raises the question of whether we can use
the same method to select remnant radio galaxies. Part of the difficulty here
is ascertaining a definition for active and remnant sources based on the
fitted remnant ratio parameter; logically, active sources would only
correspond to $R_{\mathrm{rem}}=0$ (i.e. zero off-time), however in practice
it would be impossible to distinguish a very recently switched-off source from
an active source.

By creating mock radio source catalogues with known remnant ratios, we can use
the distribution in the fitted remnant ratio to quantify the probability that:
1) an active radio source is fit with a non-zero remnant ratio; and 2) an
“unambiguous” (i.e. $R_{\mathrm{rem}}\geq 0.1$) remnant radio source is fit
with $R_{\mathrm{rem}}=0$. To do this, we follow the method outlined in Sect.
5.1 to create three mock catalogues per redshift. Each catalogue contains 5000
sources, and the remnant ratios are fixed at $R_{\mathrm{rem}}=0,0.1,0.2$,
respectively. For each redshift and source age bin, we characterise the 68th
and 95th percentiles of the fitted remnant ratio distribution; these serve as
upper bounds for the $R_{\mathrm{rem}}=0$ catalogue, and lower bounds for the
$R_{\mathrm{rem}}\geq 0.1$ catalogues. An example of this is demonstrated in
Fig. 9 for the highest source age bin at $z=0.1$, and we summarise our results
in Tab. 8 for all the redshift and source age bins.

We find that the recovery rate of active radio galaxies is high; in $95\%$ of
the cases, the fitted remnant ratio is $R_{\mathrm{rem}}\leq 0.1$ (the actual
value being 0%), which is true across all redshift and source age bins (column
4 of Tab. 8). We conclude similar results for the $R_{\mathrm{rem}}=0.1$
remnant radio galaxies; in $95\%$ of the cases, the fitted remnant ratio is
$R_{\mathrm{rem}}>0.01$ (column 6 of Tab. 8), meaning that active sources are
not contaminated by $R_{\mathrm{rem}}=0.1$ remnants, at this level. For
example, if examining a sub-catalogue of $z=0.1$ radio source with fitted ages
of 133 to 316 Myrs, we can confidently (at $2\sigma$ level) say that those
with remnant ratio $R_{\mathrm{rem}}<0.03$ will not have significant
contamination from “unambiguous”-remnants and can therefore be classified as
active sources. Conversely, those objects with fitted remnant ratios
$R_{\mathrm{rem}}>0.05$ will not have significant contamination from active
sources and can be classified as remnants. We examine the same statistics for
the $R_{\mathrm{rem}}=0.2$ sample, and reassuringly find no contamination of
active sources from this population.

##  6 CONCLUSIONS

This paper investigates a new technique to constrain the energetics of remnant
radio galaxies based on the backflow of their lobe plasma. We employ
observations of the remnant radio galaxy J2253-34 together with mock radio
source populations, in order to verify the method and examine its applications
and limitations. The approach taken to explore this technique, and the
insights gathered, are summarised below.

First, we use established spectral ageing techniques to derive a precise
measurement of the off-time in J2253-34; this is needed to independently
validate the results of our new method. To do this, we collect new radio
observations using the MeerKAT, ATCA and uGMRT telescopes (Sect. 2). The
injection index is tightly constrained by jointly fitting the JP spectra
arising from narrow regions across the lobes (Sect. 3.4.1). We then make
precise measurements of the break frequency and remnant ratio by fitting the
TCI model to the integrated radio lobe spectrum (Sect. 3.4.2); here, the
previously-fitted injection index is used as a model input. Our key results
are summarised as follows:

  * •

We develop a python implementation of synchrotron spectral ageing models, as
well as algorithms to fit these models to observed radio data. We release this
code, the synchrofit library, to the broader astrophysics community for its
relevance to modelling the synchrotron spectra from active galaxies and
supernovae remnants (Sect. 3.2).

  * •

To place the tightest constraints on the break frequency and remnant ratio,
the injection index should be independently fit and used as an input on the CI
model; this can be done by considering individual JP spectra across the lobes
(cf. Shulevski et al., 2017).

  * •

The shape of the Tribble spectral ageing models is (at most) weakly dependent
on the average magnetic field of the underlying Maxwell-Boltzmann distribution
(Sect. 3.3). These models, which likely provide a more realistic treatment of
a turbulent magnetic field, can therefore be fit to observed radio data
without needing prior knowledge of the lobe magnetic field strength.

  * •

The remnant ratio in J2253-34 is tightly constrained as
$R_{\mathrm{rem}}=0.23\pm 0.02$.

Our new method utilises the lobe surface brightness distribution to
parameterise the backflow of plasma, and exploits the RAiSE dynamics and
emissivity model. The high-resolution hydrodynamical simulation underpinning
the RAiSE code enables the synthesis of mock surface-brightness maps, which
are used to measure predicted radio source attributes. These are compared with
observations so that the intrinsic source properties can be obtained via a
parameter inversion. We apply this method to constrain the energetics of
J2253-34, and obtain the following results:

  * •

We propose two new attributes to quantify the spectral ageing of lobe plasma
as it flows away from the hotspot towards the core. The “extent 1” attribute
considers how far back the lobe visibly sweeps towards the core, and the
“extent 2” attribute considers the surface brightness of the lobe at the core,
relative to the peak surface brightness.

  * •

We implement a RAiSE-based parameter inversion using these two new attributes,
in addition to the lobe radio luminosity, size, and axis ratio, to constrain a
remnant ratio of $R_{\mathrm{rem}}=0.26\pm 0.02$ for J2253-34. Importantly,
this is consistent at the $1\sigma$ level with that constrained by our
spectral fitting, verifying our approach to constrain the dynamics of the
remnant phase based on the surface brightness distribution of the lobe
backflow.

  * •

By comparing different sets of attributes, we find that the break frequency is
the dominates the “extent” attributes fitting of the energetics; particularly
the equipartition factor and remnant ratio. By excluding the break frequency
and extent attributes from the fitting, we find that the spectral index and
curvature, alone, do not offer enough constraining power to uniquely fit the
source energetics.

  * •

For J2253-34, we accurately measure a two-sided jet kinetic power of
$Q=3.98^{+0.36}_{-0.40}\times 10^{38}$ W, a lobe magnetic field strength of
$B=0.20\pm 0.02$ nT, and importantly an active and remnant timescale of
$t_{\mathrm{on}}=37.2\pm 4.3$ Myr and $t_{\mathrm{rem}}=13.1\pm 1.5$ Myr.

We then used synthetic radio lobes, simulated for a typical 150 MHz LOFAR-like
and 1.4 GHz MeerKAT-like observation, to further examine the limitations of
our proposed method with respect to intrinsic radio source parameters. Our key
results are summarised as follows:

  * •

The jet power, source age, equipartition factor, and remnant ratio are
consistently recovered, noting that the former two are constrained with high
precision. This result builds on the work conducted by Turner et al. (2018b)
for active sources only, who similarly find that the jet power, source age,
and equipartition factor are constrained from the source size, radio
luminosity, and break frequency.

  * •

Our results show that that the parameters are more accurately fit for radio
sources with: higher redshifts, due to increased strength of IC losses;
higher-powered jets, due to increase in magnetic field strength; and for a
greater duration of the remnant phase. We therefore find that the constraining
power of our method increases when the magnitude of the energy losses,
sustained during the remnant phase, are greater.

  * •

Remnant ratios fitted by this method can be trusted if the measured source age
is large, and conversely, is less accurate for younger sources. To interpret
the on/off-times fit by this method, one should first consider the age of the
source to determine if the fitted remnant ratio is estimated reliably.

  * •

We find that our method is capable of selecting active and remnant radio
galaxy candidates. At low redshift ($z\sim 0.1$), we find that sufficiently
aged sources ($\tau>50$ Myr) fitted with $R_{\mathrm{rem}}\leq 0.1$ are active
radio galaxies within a confidence of $95\%$; this extends to sources as young
as $\tau\sim 10$ Myr at higher redshifts ($z\gtrsim 0.5$). Similarly, we find
that sources fit with $R_{\mathrm{rem}}\geq 0.1$ are true remnant radio
galaxies within a confidence of $95\%$; these classifications are more robust
at higher redshift.

The outcomes of this work show promising applications to studying the
energetic impact of radio-loud AGNs on their environments, in particular,
through studying remnant radio galaxies. However, although this method is
clearly successful for remnants with $R_{\mathrm{rem}}\lesssim 0.23$, for much
older remnants the model assumptions in the version of RAiSE used in this work
may no longer capture all the necessary physics, for example buoyancy and
Rayleigh-Taylor instability; capturing these physical processes will be an
important step necessary for modelling aged remnants, e.g. Blob 1 with
$R_{\mathrm{rem}}\sim 0.8$. Other complexities, such as dynamic or non
cluster-centered environments, can lead to asymmetries in the radio lobe
morphologies; e.g. in jet lengths, as well as lobe luminosities.
Hydrodynamical simulations of the jet and lobe evolution in such environments
already exist (see; Yates-Jones et al., 2021), which can incorporate the
calculation of the resolved adiabatic and radiative loss processes following
Yates-Jones et al. (2022). While the direct outputs of such hydrodynamical
simulations can be implemented in the grid-based approach of this work, their
synthesis is computationally expensive and can realistically only coarsely
sample any given parameter space. The RAiSE framework can naturally consider
asymmetries by allowing its governing differential equations for the source
expansion to be applied on arbitrarily small solid angle elements; in this
manner, future work can utilise modified versions of this code that include
relevant physical processes to efficiently simulate millions of mock sources
densely throughout the intrinsic parameter space.

## ACKNOWLEDGEMENTS

BQ acknowledges a Doctoral Scholarship and an Australian Government 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.
We acknowledge the Pawsey Supercomputing Centre which is supported by the
Western Australian and Australian Governments. The Australia Telescope Compact
Array is part of the Australia Telescope National Facility which is funded by
the Australian Government for operation as a National Facility managed by
CSIRO. We thank the staff of the GMRT that made these observations possible.
GMRT is run by the National Centre for Radio Astrophysics of the Tata
Institute of Fundamental Research. We additionally would like to thank
Aleksandar Shulevski for his helpful suggestions that have improved our
manuscript.

## DATA AVAILABILITY

The raw data products underlying Sect. 2 are available on the following data
archives: the Australia Telescope Online Archive (ATOA) at
https://atoa.atnf.csiro.au/ (Sect. 1); the South African Radio Astronomy
Observatory (SARAO) Web Archive at https://archive.sarao.ac.za/ (Sect. 2.2);
and the GMRT Online Archive (GOA) at
https://naps.ncra.tifr.res.in/goa/data/search (Sect. 2.3). Processed data
products underlying this article will be shared on reasonable request to the
authors. Access to the code-bases used throughout this work are disclosed
within the paper.

Table 9: Summary of the integrated radio spectrum compiled for J2253-34. Telescope | Type | Frequency | Bandwidth | Integrated flux density
---|---|---|---|---
(MHz) | (MHz) | (mJy)
MWA Phase I | Survey (GLEAM) | 119.04 | 30.72 | $238\pm 35$
154.88 | $194\pm 35$
uGMRT | Pointed observations | 417 | 250 | $114.00\pm 8.88$
682 | $79.8\pm 5.2$
ASKAP | Survey (EMU early-science) | 887 | 288 | $64.1\pm 3.2$
MeerKAT | Pointed observations | 911 | 42 | $63.17\pm 1.89$
1022 | 120 | $57.44\pm 1.85$
1422 | 153 | $41.28\pm 1.34$
1653 | 30 | $35.28\pm 1.08$
VLA | Survey (NVSS) | 1400 | 50 | $42.4\pm 3.05$
ATCA | Pointed observations | 2000 | 500 | $29.58\pm 1.77$
2868 | $19.60\pm 1.17$
4782 | $10.50\pm 0.42$
5243 | $9.38\pm 0.37$
5746 | $7.93\pm 0.32$
6216 | $6.84\pm 0.27$
8284 | $4.13\pm 0.19$
8726 | $3.84\pm 0.16$
9253 | $3.41\pm 0.15$
9716 | $3.08\pm 0.14$

## References

  * Akaike (1974) Akaike H., 1974, IEEE Transactions on Automatic Control, 19, 716
  * Alexander & Hickox (2012) Alexander D. M., Hickox R. C., 2012, New A Rev., 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
  * Binney et al. (2007) Binney J., Alouani Bibi F., Omma H., 2007, MNRAS, 377, 142
  * 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
  * Briggs (1995) Briggs D. S., 1995, in American Astronomical Society Meeting Abstracts. p. 112.02
  * Churazov et al. (2001) Churazov E., Brüggen M., Kaiser C. R., Böhringer H., Forman W., 2001, ApJ, 554, 261
  * 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
  * Cordey (1987) Cordey R. A., 1987, MNRAS, 227, 695
  * Croton et al. (2006) Croton D. J., et al., 2006, MNRAS, 365, 11
  * Duchesne & Johnston-Hollitt (2019) Duchesne S. W., Johnston-Hollitt M., 2019, PASA, 36, e016
  * 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
  * Fanaroff & Riley (1974) Fanaroff B. L., Riley J. M., 1974, MNRAS, 167, 31P
  * Frater et al. (1992) Frater R. H., Brooks J. W., Whiteoak J. B., 1992, Journal of Electrical and Electronics Engineering Australia, 12, 103
  * Godfrey et al. (2017) Godfrey L. E. H., Morganti R., Brienza M., 2017, MNRAS, 471, 891
  * Gültekin et al. (2009) Gültekin K., et al., 2009, ApJ, 698, 198
  * Hardcastle (2013) Hardcastle M. J., 2013, MNRAS, 433, 3364
  * Häring & Rix (2004) Häring N., Rix H.-W., 2004, ApJ, 604, L89
  * Harwood et al. (2013) Harwood J. J., Hardcastle M. J., Croston J. H., Goodger J. L., 2013, MNRAS, 435, 3353
  * Harwood et al. (2015) Harwood J. J., Hardcastle M. J., Croston J. H., 2015, MNRAS, 454, 3403
  * 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, PASA, 36, e048
  * Huynh et al. (2015) Huynh M. T., Bell M. E., Hopkins A. M., Norris R. P., Seymour N., 2015, MNRAS, 454, 952
  * Huynh et al. (2020) Huynh M. T., Seymour N., Norris R. P., Galvin T., 2020, MNRAS, 491, 3395
  * Jaffe & Perola (1973) Jaffe W. J., Perola G. C., 1973, A&A, 26, 423
  * Jarvis et al. (2016) Jarvis M., et al., 2016, in MeerKAT Science: On the Pathway to the SKA. p. 6 (arXiv:1709.01901)
  * Johnston et al. (2007) Johnston S., et al., 2007, PASA, 24, 174
  * Jonas & MeerKAT Team (2016) Jonas J., MeerKAT Team 2016, in MeerKAT Science: On the Pathway to the SKA. p. 1
  * 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 & Alexander (1997) Kaiser C. R., Alexander P., 1997, MNRAS, 286, 215
  * Kaiser et al. (1997) Kaiser C. R., Dennett-Thorpe J., Alexander P., 1997, MNRAS, 292, 723
  * Kale & Ishwara-Chandra (2021) Kale R., Ishwara-Chandra C. H., 2021, Experimental Astronomy, 51, 95
  * Kapinska et al. (2015) Kapinska A. D., Hardcastle M., Jackson C., An T., Baan W., Jarvis M., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14). p. 173 (arXiv:1412.5884)
  * Kardashev (1962) Kardashev N. S., 1962, Soviet Ast., 6, 317
  * Komatsu et al. (2011) Komatsu E., et al., 2011, ApJS, 192, 18
  * Komissarov & Gubanov (1994) Komissarov S. S., Gubanov A. G., 1994, A&A, 285, 27
  * Lacy et al. (2020) Lacy M., et al., 2020, PASP, 132, 035001
  * Longair (2011) Longair M. S., 2011, High Energy Astrophysics
  * Magorrian et al. (1998) Magorrian J., et al., 1998, AJ, 115, 2285
  * 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
  * McMullin et al. (2007) McMullin J. P., Waters B., Schiebel D., Young W., Golap K., 2007, in Shaw R. A., Hill F., Bell D. J., eds, Astronomical Society of the Pacific Conference Series Vol. 376, Astronomical Data Analysis Software and Systems XVI. p. 127
  * McNamara & Nulsen (2007) McNamara B. R., Nulsen P. E. J., 2007, ARA&A, 45, 117
  * Miley (1980) Miley G., 1980, ARA&A, 18, 165
  * Morganti (2017) Morganti R., 2017, Nature Astronomy, 1, 596
  * Murgia et al. (2011) Murgia M., et al., 2011, A&A, 526, A148
  * Nagai et al. (2006) Nagai H., Inoue M., Asada K., Kameno S., Doi A., 2006, ApJ, 648, 148
  * Norris (2011) Norris R. P., 2011, Journal of Astrophysics and Astronomy, 32, 599
  * Offringa & Smirnov (2017) Offringa A. R., Smirnov O., 2017, MNRAS, 471, 301
  * Offringa et al. (2012) Offringa A. R., van de Gronde J. J., Roerdink J. B. T. M., 2012, A&A, 539, A95
  * Offringa et al. (2014) Offringa A. R., et al., 2014, WSClean: Widefield interferometric imager (ascl:1408.023)
  * Oosterloo et al. (2010) Oosterloo T., Verheijen M., van Cappellen W., 2010, in van Leeuwen J., ed., ISKAF2010 Science Meeting. p. 43 (arXiv:1007.5141)
  * Pacholczyk (1970) Pacholczyk A. G., 1970, in Series of Books in Astronomy and Astrophysics.
  * Parma et al. (2007) Parma P., Murgia M., de Ruiter H. R., Fanti R., Mack K. H., Govoni F., 2007, A&A, 470, 875
  * Quici et al. (2021) Quici B., et al., 2021, PASA, 38, e008
  * Reynolds (1994) Reynolds J. E., 1994, ATNF Internal Memo AT/39.3/040, A Revised Flux Scale for the AT Compact Array. Australia Telescope National Facility
  * Rybicki & Lightman (1979) Rybicki G. B., Lightman A. P., 1979, Radiative processes in astrophysics
  * Sault et al. (1995) Sault R. J., Teuben P. J., Wright M. C. H., 1995, in Shaw R. A., Payne H. E., Hayes J. J. E., eds, Astronomical Society of the Pacific Conference Series Vol. 77, Astronomical Data Analysis Software and Systems IV. p. 433 (arXiv:astro-ph/0612759)
  * Schoenmakers et al. (2000) Schoenmakers A. P., de Bruyn A. G., Röttgering H. J. A., van der Laan H., Kaiser C. R., 2000, MNRAS, 315, 371
  * 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
  * Shulevski et al. (2012) Shulevski A., Morganti R., Oosterloo T., Struve C., 2012, A&A, 545, A91
  * Shulevski et al. (2017) Shulevski A., et al., 2017, A&A, 600, A65
  * Tadhunter (2008) Tadhunter C., 2008, Mem. Soc. Astron. Italiana, 79, 1205
  * Tribble (1991) Tribble P. C., 1991, MNRAS, 253, 147
  * 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 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
  * Vikhlinin et al. (2006) Vikhlinin A., Kravtsov A., Forman W., Jones C., Markevitch M., Murray S. S., Van Speybroeck L., 2006, ApJ, 640, 691
  * Worrall et al. (2012) Worrall D. M., Birkinshaw M., Young A. J., Momtahan K., Fosbury R. A. E., Morganti R., Tadhunter C. N., Verdoes Kleijn G., 2012, MNRAS, 424, 1346
  * Yates-Jones et al. (2021) Yates-Jones P. M., Shabala S. S., Krause M. G. H., 2021, MNRAS, 508, 5239
  * Yates-Jones et al. (2022) Yates-Jones P. M., Turner R. J., Shabala S. S., Krause M. G. H., 2022, MNRAS, 511, 5225
  * van Haarlem et al. (2013) van Haarlem M. P., et al., 2013, A&A, 556, A2