1. Introduction
Pairs of radio-emitting jets launched via accretion processes in radio-loud active galactic nuclei (RLAGN) are among the most energetic objects in the Universe and are known to be an episodic phenomenon, operating on a duty cycle over the lifetime of the host galaxy. The most obvious evidence for episodic jet activity is shown by ‘double-double’ radio galaxies (DDRGS; Schoenmakers et al. Reference Schoenmakers, de Bruyn, Rottgering, van der Laan and Kaiser2000) in which two pairs of radio lobes are observed; the inner pair signalling recent jet production and the outer pair leftover from a past cycle of activity. Further direct evidence is supplied by the ripples and multiple X-ray cavities seen in the Perseus (Boehringer et al. Reference Boehringer, Voges, Fabian, Edge and Neumann1993; Fabian et al. Reference Fabian2003) and Virgo clusters (Forman et al. Reference Forman2005). Indirect evidence of a strong correlation between quiescent star formation and the presence of active RLAGN in local, massive galaxies implies prolonged radio-loud phases where the emitted energy is sufficient to prevent the cooling of halo gas (Best et al. Reference Best2005, Reference Best, Von Der Linden, Kauffmann, Heckman and Kaiser2007; Barišić et al. Reference Barišić2017; Sabater et al. Reference Sabater2019). From a theoretical perspective, episodic behaviour is expected and necessary for a self-regulating feedback loop to achieve a balance between cooling, infalling gas, and the heating, and redistribution of that gas (e.g. Kawata & Gibson Reference Kawata and Gibson2005; Shabala & Alexander Reference Shabala and Alexander2009a; Novak, Ostriker, & Ciotti Reference Novak, Ostriker and Ciotti2011; Pope, Mendel, & Shabala Reference Pope, Mendel and Shabala2012; Raouf et al. Reference Raouf, Shabala, Croton, Khosroshahi and Bernyk2017).
The lifecycle of an AGN is often broken down into four stages in the literature; the youth stage, the evolved active phase, the remnant phase, and the restarted phase. A schematic of the AGN lifecycle is shown in Figure 1. The youth stage is often identified with Compact Steep Spectrum (CSS) or Gigahertz Peaked Spectrum (GPS) sources (O’Dea Reference O’Dea1998). As they pass through the host galaxy, RLAGN jets drive gas outflows (Morganti, Tadhunter, & Oosterloo Reference Morganti, Tadhunter and Oosterloo2005; Morganti et al. Reference Morganti, Fogasy, Paragi, Oosterloo and Orienti2013; Mahony et al. Reference Mahony2015; Santoro et al. Reference Santoro, Tadhunter, Baron, Morganti and Holt2020), create turbulence, and shock-heat the ISM, affecting star formation rates (either by enhancing it as in Croft et al. Reference Croft2006; Dugan, Gaibler, & Silk Reference Dugan, Gaibler and Silk2017; Nesvadba et al. Reference Nesvadba, Bicknell, Mukherjee and Wagner2020; Zovaro et al. Reference Zovaro2020, or by suppressing it as in Nesvadba et al. Reference Nesvadba2010; Nesvadba et al. Reference Nesvadba, Boulanger, Lehnert, Guillard and Salome2011; Morganti et al. Reference Morganti, Fogasy, Paragi, Oosterloo and Orienti2013; Mandal et al. Reference Mandal2021).
A schematic of the AGN lifecycle from an active source (left), to a remnant source (middle) to restarted sources (right). If the time for the remnant lobes to fade below the detectable limit is longer than the time for the nuclear activity to reignite, then remnant lobes may be seen together with a newly restarted, young radio source (bottom right). Else, a young radio source may be seen with no indication of past activity (top right).
If the RLAGN jets evolve to larger scales,Footnote a the jet-inflated lobes can modify the thermal state of the surrounding environment, quenching otherwise rapid cooling at the centres of galaxy clusters (Churazov et al. Reference Churazov, Bruggen, Kaiser, Bohringer and Forman2001; Shabala & Alexander Reference Shabala and Alexander2009b; Mittal et al. Reference Mittal, Hudson, Reiprich and Clarke2009; Fabian Reference Fabian2012; Yang & Reynolds Reference Yang and Reynolds2016). The effects of ‘jet-mode’ feedback are now explicitly included in all modern cosmological galaxy formation models (Sijacki et al. Reference Sijacki, Springel, Di Matteo and Hernquist2007; Schaye et al. Reference Schaye2014; Weinberger et al. Reference Weinberger2016; Davé et al. Reference Davé2019; Li et al. Reference Li2021; Bird et al. Reference Bird2022; Villaescusa-Navarro et al. Reference Villaescusa-Navarro2023). Cosmological models implement jet feedback in a variety of ways ranging from isotropic energy-momentum injection to collimated beams. However, due to the hugely dynamic range that would need to be covered, all cosmological models struggle to capture important features of real jets including relativistic jet speeds and initially wide opening angles. On the other hand, dedicated jet simulations can miss the complexity of cosmological environments. A few works have brought the two together including Mendygral, Jones, & Dolag (Reference Mendygral, Jones and Dolag2012) and the CosmoDRAGoN simulation suite of Yates-Jones et al. (Reference Yates-Jones2023).
Following an active episode, the remnant stage (the focus of this work) will begin where the nuclear activity ceases or substantially reduces (Sabater et al. Reference Sabater2019), the jet outflow stops, and the remnant lobes are left to fade and rise buoyantly out of the central gravitational well (Churazov et al. Reference Churazov, Bruggen, Kaiser, Bohringer and Forman2001; Reynolds et al. Reference Reynolds, McKernan, Fabian, Stone and Vernaleo2005). As they rise, they entrain and uplift large quantities of low entropy gas to large radii (Chen, Heinz, & Enßlin Reference Chen, Heinz and Enßlin2019). Finally, a restarted phase may occur where the nuclear activity reignites and a new jet pair may be observed (Brocksopp et al. Reference Brocksopp, Kaiser, Schoenmakers and de Bruyn2011; Orrù et al. Reference Orrù2015; Nandi et al. Reference Nandi2019; Mahatma Reference Mahatma2023).
Compared to the active phase of RLAGN evolution, the remnant phase is considerably less-well studied. Only a handful of studies have presented small observed samples of remnant candidates identified based on the assumptions about the morphological and spectral properties of aged radio plasma (e.g. Murgia et al. Reference Murgia2011; Mahatma et al. Reference Mahatma2018; Jurlin et al. Reference Jurlin2020; Jurlin et al. Reference Jurlin2021; Quici et al. Reference Quici2021; Dutta et al. Reference Dutta2023). The numerical modelling that has concentrated on the large-scale evolution of the remnant phase has been mostly focused on the energy transfer between the remnant lobe and surrounding medium for powerful (Q
$_{\mathrm{j}} = 10^{38-49}$
W), large-scale (
$\gt100$
kpc), and, often, FR-II-like sources in smooth, radially symmetric environments (Zanni et al. Reference Zanni2005; Perucho, Quilis, & Martí Reference Perucho, Quilis and Mart2011; Perucho et al. Reference Perucho, Mart, Quilis and Ricciardelli2014; Walg et al. Reference Walg, Achterberg, Markoff, Keppens and Porth2014; English, Hardcastle, & Krause Reference English, Hardcastle and Krause2019; Chen et al. Reference Chen, Heinz and Enßlin2019). Similarly, analytical models have focused on describing the transition of powerful, strongly overpressured active sources (e.g. Kaiser & Cotter Reference Kaiser and Cotter2002; Hardcastle Reference Hardcastle2018). These studies describe a system where, after the jets switch off, the radio lobes continue to rise through the ambient atmosphere, potentially driven by forces beyond simple buoyancy (e.g. Perucho et al. Reference Perucho, Mart, Quilis and Ricciardelli2014; English et al. Reference English, Hardcastle and Krause2019). During the transition from active to remnant source, a significant portion of the lobe energy may remain contained within the lobes at the end of the active phase (English et al. Reference English, Hardcastle and Krause2019). Various mechanisms for energy transfer have been proposed. Entropy studies of the surrounding environment by Zanni et al. (Reference Zanni2005) highlight bow shock heating as the primary driver. In contrast, Chen et al. (Reference Chen, Heinz and Enßlin2019) argue that the thermal exchange between uplifted, low-entropy gas and the ambient medium plays a significant role in heating cool-core systems. The overall paucity of remnant studies is unfortunate, as these systems hold key information about the late-stage evolution of RLAGN, the efficiency of energy injection into the surrounding medium, and the timescales over which feedback persists (Turner Reference Turner2018; Quici et al. Reference Quici, Turner, Seymour and Hurley-Walker2025).
In a previous work, (Yates-Jones et al. Reference Yates-Jones2023, hereafter Paper I) we presented the first results from the CosmoDRAGoN project; a large suite of AGN jet simulations carried out in environments derived from cosmological simulations. In that paper, the dynamic properties and spatially resolved spectral properties of active sources were primarily considered for simulations spanning both low and high jet kinetic powers corresponding to FR-I and FR-II morphologies. The evolution of remnants was only briefly considered.
The present paper extends the results of Paper I by looking more closely at the transition from the active to remnant phase. We choose to use the CosmoDRAGoN suite of simulations, as natural inhomogeneities in a dynamic environment are likely to become important for the spatial distribution of remnant plasma and the timescales over which instabilities develop. Using these more complex environments, we explore the mechanisms driving lobe evolution after the active phase ends and compare the results from our numerical simulations to analytical expectations. In this paper, our results are based on data obtained from the raw hydrodynamic outputs. A study of the spectral properties of the simulations is deferred to the companion paper (Stewart et al. in preparation). The current paper is organised as follows: Section 2 describes the computational setup and input parameters for our simulation suite. In Section 3 we present our findings, and in Section 4, we provide a discussion of the results including comparisons to analytic modelling. Finally, in Section 5 we summarise the work and present our conclusions.
2. Simulation setup and parameters
To compile our simulation suite, we take five existing active simulations (i.e. remnant progenitors) from the CosmoDRAGoN simulation catalogue presented in Paper I and continue their evolution in the remnant phase. In Section 2.1 we outline the numerical techniques used for all simulations. We then describe the method for environment selection and jet injection in Sections 2.2 and 2.3, respectively. Because this work models the transition from active to remnant radio sources, we explain the process to terminate the jet fluid and evolve the remnant in Section 2.4. Finally, we describe our choice of remnant progenitor properties in Section 2.5.
2.1. Numerical methods
All active and remnant simulations used in this work are run using a modified version of the publicly-available pluto code version 4.3.Footnote
b
The numerical techniques employed follow those of Yates-Jones, Shabala, & Krause (Reference Yates-Jones, Shabala and Krause2021) and Paper I. The fluid variables are related by the Taub-Mathews equation of state (Taub Reference Taub1948; Mathews Reference Mathews1971; Mignone & McKinney Reference Mignone and McKinney2007) and we solve the Eulerian equations of relativistic hydrodynamics using the hllc Riemann solver. We use a linear reconstruction scheme and a Courant-Friedrichs-Lewy (CFL) number of 0.3. A second-order Runge-Kutta integration progresses the simulation in time and a minmod flux limiter is used to prevent spurious oscillations in the presence of shocks and discontinuities. We adopt the standard
$\Lambda$
CDM cosmology, which is consistent with the simulation catalogue that our environments are derived from (see Section 2.2 below). The parameters are taken from the Planck mission and are
$\Omega_{\textrm{M}} = 0.307$
,
$\Omega_{\textrm{B}} = 0.048$
,
$\Omega_{{\Lambda}} = 0.693$
(Planck Collaboration et al. 2016).
All simulations are run on a static, three-dimensional Cartesian domain (X, Y, Z) with maximum physical dimensions: 400 kpc
$\times 400$
kpc
$\times 400$
kpc (−200 kpc to +200 kpc in all dimensions). The grid extents for a given simulation are chosen based on the expected size of the simulation (some simulations in this work do not require the full 400 kpc grid) but the resolution remains consistent across all runs.
All three dimensions are separated into five resolution patches such that the finest resolution is only used where it is needed (i.e. near the jet injection region). The inner (
$-2.5$
kpc,
$2.5$
kpc) patch is uniformly covered by 100 cells giving a resolution of 0.05 kpc/cell. In the intermediate (
$\pm 2.5$
kpc,
$\pm 10$
kpc) region, a geometrically stretched grid is applied over 100 cells with a stretching ratio of
$1.0076$
. This gives a resolution of 0.08 kpc/cell at 10kpc. In the outer (
$\pm 10$
kpc,
$\pm 200$
kpc) patch, a stretching ratio of
$1.0084$
is used with a resolution of 0.84 kpc/cell at 100 kpc. These stretching ratios result in a total grid cell count of 960 cells in each dimension. Following Paper I, we impose outflow boundary conditions at all grid edges to avoid amplifying any initially small inhomogeneities present in our cosmological environments (see Section 2.2).
The simulations of the active phase were run on the Gadi facility provided by the National Computational Infrastructure, Australia. The follow-up remnant simulations were carried out using the kunanyi facility provided by the Tasmanian Partnership for Advanced Computing. For our most computationally intensive model, the total computational time was approximately 790 000 CPU hours. The primary data products output from the simulation are the hydrodynamic quantities of density, pressure, velocity (in three dimensions), and a passive jet tracer value for each cell at each snapshot.
2.2. Environment implementation
The environment in which a radio galaxy evolves can have a significant impact on the morphology, dynamics, and observable properties such as the size, shape, brightness and spectral features of the radio emission (Kaiser & Alexander Reference Kaiser and Alexander1997; Rodman et al. Reference Rodman2018; Yates-Jones et al. Reference Yates-Jones, Shabala and Krause2021). The properties of this simulated environment are a particularly important consideration of the current study given that, in the absence of a driving momentum flux, remnant plasma will be more susceptible to the properties and inhomogeneities of the surrounding medium than an active source.
Some remnants have been associated with both cluster and group environments, while other remnants are considered to be isolated (Murgia et al. Reference Murgia2011; Jurlin et al. Reference Jurlin2021). In our study, we choose to analyse simulations run in complex cluster and group-like environments. These complex environments are extracted from the pressure, temperature, and density profiles from the the three hundred project catalogue of re-simulated galaxy clusters with full-physics hydrodynamics (Cui et al. Reference Cui2018). The process for doing this is outlined in-depth in Paper I, but we provide a brief overview here.
First, suitable halo candidates are selected from the the three hundred project catalogue. These are selected based on size, and overall stability, ensuring that there is no evidence of merger events at the epoch of interest. We use data from the
$z=0$
snapshots are used throughout this paper.
Second, the pressure, density, momentum density, and force density are interpolated onto a three-dimensional Cartesian grid with a
$1\,\mathrm{kpc}$
/cell resolution such that total mass is conserved. The standard smooth particle hydrodynamics (SPH) ‘scatter’ formalism, implemented in sphtool
Footnote
c
, is used to smooth the environment quantities, such that the smoothed quantity
$A_s$
is given as a function of position
r
and the SPH properties of i particles,
\begin{equation} A_s(\mathbf{r}) \approx \sum_i m_i \frac{A_i}{\unicode{x03C1}_i}W(|\mathbf{r}-\mathbf{r}_{\mathbf{i}}|, h_i) \end{equation}
with particle mass, density, and smoothing length
$m_i$
,
$\unicode{x03C1}_i$
, and
$h_i$
, respectively. We use the
$M_4$
kernel (Monaghan & Lattanzio Reference Monaghan and Lattanzio1985) for the smoothing function,
$W(\mathbf{r}, h)$
, as detailed further in Paper I. The velocity and acceleration fields are, respectively, derived from the momentum and force fields. As pluto version 4.3 does not support self-gravity, gadget-2 is used to calculate the gravitational acceleration necessary to generate gravitationally stable initial conditions. The gravitational acceleration is then interpolated using the same method as for the hydrodynamic quantities and initialised on the grid as a static field. This field cannot be evolved throughout the pluto simulations, resulting in a finite time before the gas density and gravitational potential fields dissociate from each other. However, stability tests conducted in Paper I show the environments used in this work to be stable for at least 300 Myr, which is longer than the timescales for which we simulate our remnant radio galaxies.
In the present work, we use remnant progenitor simulations run in halos 0003 and 0031 from cluster 002 in the three hundred project catalogue as the respective cluster- and group-like environments. These environments are chosen based on their long-term stability. A list of parameters for these halos is provided in Table 1 and in Figure 2 we show the ambient density and pressure profiles.
Cosmological environment values for central density, central pressure, halo mass, virial radius, and average core radius.
Environment density (left) and pressure (right) profiles are shown for the cosmological cluster (orange) and group (green) environments. Respectively, the dashed, dot-dashed and dotted lines show the radial evolution along the z-, x- and y- axes. The solid lines indicate the symmetric, radially averaged fit to the cosmological environments. For comparison, the range of environment profiles considered by the related work of English et al. (Reference English, Hardcastle and Krause2019) is shown by the shaded region.
2.3. Jet injection
The remnant progenitors chosen from the CosmoDRAGoN simulation suite are injected as a pair of bipolar, conical outflows with half-opening angle
$\unicode{x03B8}_{\text{j}}$
. They are injected through a spherical, internal boundary condition of radius
$r_{\textrm{o}}=0.75$
kpc centred at the origin. Within
$r \lt r_{\textrm{o}}$
, the pressure and density variables are continuously overwritten with the input jet density,
$\unicode{x03C1}_{\mathrm{j}}$
, and pressure
$\mathrm{p}_{\mathrm{j}}$
which are calculated following
\begin{equation}\unicode{x03C1} (r \lt r_0) = 2\unicode{x03C1}_{\mathrm{j}}\left[1+\left(\frac{r}{r_0}\right)^2\right]^{-1} \end{equation}
\begin{equation}P (r \lt r_0) = 2^\Gamma P_{\text{j}} \left( \frac{\unicode{x03C1}(r)}{2\unicode{x03C1}(r_{0})}\right)^\Gamma, \end{equation}
where r is the spherical radius from the grid origin and
$\Gamma$
is the ideal gas adiabatic index calculated according to the Taub-Mathews equation of state. Within the jet injection cone,
$\unicode{x03B8} \leq \unicode{x03B8}_j$
, the fluid velocity is set to the jet velocity,
$v_{\mathrm{j}}$
. For a relativistic regime, the one-sided jet kinetic power is related to the initial jet pressure, density, velocity, and cross-sectional area of the jet head at injection,
$A_{\mathrm{j}}$
by
\begin{equation}Q_j = \left[\frac{\Gamma}{\Gamma - 1} P_{\mathrm{j}} \unicode{x03B3}^2 + \unicode{x03B3}(\unicode{x03B3}-1)\unicode{x03C1}_{\mathrm{j}} c^2 \right]v_{\mathrm{j}} A_{\mathrm{j}}, \end{equation}
where c is the speed of light in a vacuum and
$\unicode{x03B3} = 1/\sqrt{1-v_{\mathrm{j}}^2/c^2}$
is the bulk flow Lorentz factor. For non-relativistic, strongly supersonic flows, equation (4) reduces to the flux of kinetic energy density along the jet,
\begin{equation}Q = \frac{1}{2}\unicode{x03C1}_{\mathrm{j}} v_{\mathrm{j}}^3A_{{j} }. \end{equation}
A passive tracer fluid is injected with the jet material, taking an initial value of 1 within the injection region and 0 elsewhere.
2.4. Jet termination
In this study, we focus on the evolution of simulated active radio jets into the remnant phase. To do this, we terminate the constant flow of injected material through the injection region and continue the simulation, following the evolution of the remnant lobes. When the jet flow is stopped, the fluid variables in the spherical injection region are no longer overwritten and assume similar values to their surrounding cells. To mitigate any adverse numerical affects after the jet is switched off, we include a ‘ramp-down’ of the jet velocity over a time period of 0.05 Myr.
The timescales over which active sources transition to remnants are unclear. Optical and radio investigations suggest that the active phase can occur over long periods of 10–100 Myr (Parma et al. Reference Parma2007; Murgia et al. Reference Murgia2011; Hardcastle et al. Reference Hardcastle2019; Shabala et al. Reference Shabala2020) as well as shorter periods of 0.01–0.1 Myr (Murgia Reference Murgia2003; Polatidis & Conway Reference Polatidis and Conway2003). Once the jet flow is terminated, the remnant phase has been estimated to last anywhere from less than half of the active phase (Parma et al. Reference Parma2007; Orrù et al. Reference Orrù2010; Murgia et al. Reference Murgia2011), to timescales comparable to Shulevski et al. (Reference Shulevski2015), or even longer than the active phase (Brienza et al. Reference Brienza2016)Footnote
d
but recent modelling suggests that the fading time of remnant sources is rapid (Yates, Shabala, & Krause Reference Yates, Shabala and Krause2018; English et al. Reference English, Hardcastle and Krause2019; Shabala et al. Reference Shabala2020). To explore a parameter space of varying active lifetimes and source sizes, each progenitor model (described below in Section 2.5) is terminated at three points: when the total source length has reached 20, 60, and 180 kpc. The total length of each lobe is approximately half the total source size, as all simulations show negligible asymmetries about the x-axis in the active phase. Terminating each model at 20, 60, and 180 kpc results in active lifetimes,
$t_{\textrm{on}}$
, that span 0.5–100 Myr.
2.5. Simulation parameters
Our goal is to include a range of large-scale progenitor morphologies for evolving our remnants. Our parameter space, including jet-injection properties and environment, is consistent with this aim. FR-I morphologies are thought to be the result of a combination of non-relativistic or mildly relativistic velocities (0.01-0.5c),Footnote
e
low kinetic jet powers (
$10^{36}$
–
$10^{37}$
W), and wide half-opening angles (
$\gt 24^{\textrm{o}}$
) while FR-II-like morphologies are often characterised by relativistic injection velocities, narrow half-opening angles (
$\lt24^{\textrm{o}}$
) and high kinetic powers (
$\gt10^{37}$
W) (Krause et al. Reference Krause, Alexander, Riley and Hopton2012; Massaglia et al. Reference Massaglia, Bodo, Rossi, Capetti and Mignone2016).
We explore a parameter space of low (
$Q_{\mathrm{j}} = 10^{36}$
W) and high (
$Q_{\mathrm{j}} = 10^{38}$
W) jet kinetic powers, relativistic (0.98c) and sub-relativistic (0.01c) jet injection velocities, narrow (
$\unicode{x03B8}_{\mathrm{j}} = 7.5^{\textrm{o}}$
) and wide (
$\unicode{x03B8}_{\mathrm{j}} = 25^{\textrm{o}}$
) half-opening angles across a range of active times in both group and cluster environments. The parameters for all simulation runs are shown in Table 2. Each simulation is assigned an identification code of the form Qaa-vbb-acc-Ddd where Qaa denotes the jet power (in log W), vbb is the injected velocity (as a fraction of c; i.e. v98 is 0.98c), acc is the half-opening angle (in degrees), D denotes the type of environment (C for cluster or G for group), and dd is the total source length at the start of the remnant phase (in kpc). We consider a total of 15 simulations, nine of which are run in the cluster environment, and six in the group.
Parameters for the simulation runs.
$Q_j$
is the one-sided kinetic jet power,
$v_j$
is the initial jet velocity,
$\unicode{x03B8}_j$
is the half-opening angle,
$t_{\textrm{on}}$
is the time at which the remnant phase commences, and
$t_{\textrm{on}} + t_{\textrm{off}}$
is the total simulation time. The run code of each model is given in the last column.
3. Results
We have performed a series of hydrodynamic simulations in realistic environments to study how active and remnant dynamics are affected by different environment and jet initial conditions. We are particularly interested in the evolving morphology of the remnant lobes, their deceleration, and transition to pressure balance with the ambient medium. In this section, we examine the key morphological features of all simulations in Section 3.1. We quantify the dynamic evolution of the radio lobes in Section 3.2 before tracking the evolution of the shocked material surrounding the lobes in Section 3.3.
Mid-plane slices in the
$X-Z$
plane of the logarithmic density distributions for all models at the last active output. The total source age is displayed on each panel. From left to right, the columns show 20, 60, and 180 kpc switch-off simulations. The spatial scales have been adjusted across the three columns such that each simulation can be seen clearly. We have grouped the simulations such that low-power sources are shown in the top two rows and high-power sources in the bottom three rows. Group simulations are shown in the panel immediately above the equivalent cluster simulation. The blue contours outline where the passive tracer is above the threshold of
$10^{-4}$
.
3.1. General features of radio source evolution
3.1.1. Active phase
The morphology of each simulation is best understood by looking at the two-dimensional logarithmic density distributions. Mid-plane slices of this distribution for all simulations during the active phase are shown in Figure 3. Each panel in Figure 3 corresponds to the last data output before the end of the active phase,
$t_{\textrm{on}}$
. The values for
$t_{\textrm{on}}$
are given in Table 2. The lobes are identified by the regions of underdense material (with respect to the local ambient medium) surrounded by a dense shell of hot, shocked ambient material. These lobes map well to the passive tracer quantity (the blue contour) which outlines all material with a tracer value greater than 10
$^{-4}$
. The amount of tracer material in the lobe is largely insensitive to this value, provided it is not greater than approximately 10
$^{-2}$
. Similar tracer cutoffs are used in other related work (Hardcastle & Krause Reference Hardcastle and Krause2013; Yates et al. Reference Yates, Shabala and Krause2018 and Paper I). We consistently use a tracer cutoff of 10
$^{-4}$
throughout this paper.
The slow, low-power simulations are shown in the top two rows of Figure 3. All low-power simulations have similar jet kinetic power (
$10^{36}$
W) and injection velocity (
$0.01\textrm{c}$
) but differ in the environment they propagate into. The lobes of the low-power group (top row) and cluster (second row) simulations are filled with mildly (
$\lesssim$
0.5 dex) underdense, forward-flowing material, surrounded by an oval-shaped region of weakly shocked ambient gas. The magnitude of pressure in the shocked region is strongest at the lobe head and diminishes towards the equatorial regions. In the 180 kpc switch-off, cluster simulation, the weakly shocked material is nearly indistinguishable from the ambient medium (see simulation Q36-v01-a25-C180 in the far right panel of row two in Figure 3).
The key difference between the low-power group and cluster simulations is the spatial distribution of the lobe plasma; group simulations inflate wider lobes. All six low-power simulations are injected as initially overdense, conical outflows and are left to collimate via a recollimation shock. Since the central density and pressure in the group environment is an order of magnitude less than in the cluster environment, simulations run in the group are collimated about 3 times further downstream, and at a jet radius of around
$3.5$
times larger than the equivalent cluster simulation. This is consistent with analytic theory, where the collimation length scale goes as
$\unicode{x03C1}_{\mathrm{x}}^{-1/2}$
(Alexander Reference Alexander2006). At later times, the higher pressures and densities of the cluster environment act to confine the inflating lobe material.
Our high-power simulations are shown in the bottom three rows of Figure 3. All nine high-power simulations are run with an initial jet power of
$Q_{\mathrm{j}} = 10^{38}$
W and injection velocity of
$v_{\mathrm{j}} = 0.98$
c characteristic of powerful radio galaxies. One set of simulations have a narrow half-opening angle of
$7.5^{\rm o}$
(bottom row), the other two have a wider half-opening angle of
$25^{\rm o}$
, designed to produce FR-I sources (e.g. Krause et al. Reference Krause, Alexander, Riley and Hopton2012).
All high-power sources are collimated within 5 kpc of injection and are characterised on large scales by significantly underdense (3-4 dex lower than the ambient medium) lobes. These are primarily filled by backflow from the strong terminal shock at the head of the lobe. Surrounding the lobes, the region of shocked ambient material (the ‘shocked shell’) is roughly uniform in density. Such features are typical of overpressured outflows as seen in e.g. Komissarov & Falle (Reference Komissarov and Falle1998), Hardcastle & Krause (Reference Hardcastle and Krause2013).
For the wide jet opening angle simulations, the spatial distribution of the lobe plasma is similar regardless of environment (compare row three and four of Figure 3). However, wide opening angle simulations in the cluster environment (row four in Figure 3) consistently take twice as long as the group simulations to evolve to a certain size. For the narrow, high-power simulations, the lobes are thinner and the source takes only a third of the time to traverse the same distance as the wider angle simulation in the same environment.
3.1.2. Remnant phase
We allow each simulation to evolve for 100 Myr in the remnant phase, or until the shock front reaches the edge of the computational domain. Beyond this point, mass is no longer conserved on the grid, and the full interaction of the shock with the environment cannot be captured. Snapshots of our remnant sources 50 Myr post switch-off are given in Figure 4.
Mid-plane slices for all models 50 Myr after the remnant phase started. The layout of this figure is analogous to Figure 3.
In all sources, the shock front progressively separates from the lobes, becomes increasingly spherical, and weakens against the ambient medium. For some simulations (e.g. low-powered remnants in cluster environments), the shocked region has diminished such that it cannot be detected against the ambient medium 50 Myr after switch-off. For low-powered sources, the lobe separation increases and the inertia-driven material in the jet at switch-off continues along the evacuated channel, rising faster than the bulk lobe plasma. For our 20 kpc switch-off, low-power, cluster simulation (Q36-v01-a25-C20), the trailing jet material completely drills through the existing lobe and separates at
$\sim \pm 20$
kpc from the grid origin, forming four main components (see the left-hand panel in row two of Figure 4). For larger low-power sources, the trailing jet material runs out of thrust before it traverses the lobe, and the lobes do not become fragmented at any point during the simulated time.
The evolution of high-power, fast simulations is different. The lobes are dominated by backflow, and this continues to mix in a single mass close to the core before the lobes become pinched due to large Rayleigh-Taylor instabilities. The trailing jet material quickly (
$\lt 1$
Myr; less than the time between simulation outputs) reaches the lobe head, rebounds off the leading shock and mixes back in with older lobe material (i.e. it continues to do what the jet was doing before switch off). For large sources, the bulk flow continues in the direction of the declining ambient profile. Our 20 and 60 kpc switch-off simulations are confined to the inner 70–50 kpc for most of their active and remnant life. In this region, the ambient density only drops by a factor of two and the preferential direction of fluid flow is not along the jet axis. The evolution of these sources is then dependent on the magnitude of the velocities within the remnant lobe and the inhomogeneities in the environment. As shown in the right-hand panels in the bottom three rows of Figure 4, this results in a greater loss of structure at
$t_{\textrm{on}} + 50$
Myr. We note that the simulations presented in this paper do not contain magneto-hydrodynamic effects, which may suppress instabilities at the boundary between the lobe and shocked shell (Shabala & Alexander Reference Shabala and Alexander2009a). However, we expect the modelling of flows within the lobes to be representative as these flows are governed by pressure gradients, and large-scale hydrodynamic interactions.
3.2. Lobe dynamics
We now probe the dynamic evolution of the active and remnant lobes in more detail. We explore the evolution of models with low kinetic power progenitors (
$Q_{\mathrm{j}} = 10^{36}$
W) in Section 3.2.1 and those with high-power progenitors (
$Q_{\mathrm{j}} = 10^{38}$
W) in Section 3.2.2.
3.2.1. Low-power progenitors
The left-hand panel of Figure 5, shows the time evolution of the total source length along the jet propagation axis. This is measured from the maximum extent along the jet axis of the tracer fluid above
$10^{-4}$
. For
$t \lt t_{\textrm{on}}$
, a key feature of our low-power sources is that equivalent jet models in different environments grow at a similar rate and have similar total active times for a given size. This is shown by the proximity and slope of the thick blue and purple lines in Figure 5. This is not an intuitive result. Since the central density and pressure in the group environment are approximately an order of magnitude below the cluster, one would expect all models in the group to propagate faster. However, in the group simulations, the larger jet radius at collimation (discussed in Section 3.1.1) and the greater lateral expansion of the lobe plasma cause the momentum flux to dissipate over a larger area.
The time evolution of the total source length measured as the maximum extent (along z) of the tracer material above
$10^{-4}$
. The left-hand panel shows all models with low injected jet kinetic powers (
$Q_{\mathrm{j}} = 10^{36}$
W). Models with high injected jet kinetic powers (
$Q_{\mathrm{j}} = 10^{38}$
W) are shown on the right. The active and remnant phases are shown with thick and thin lines, respectively. The cross-marks indicate
$2\times t_{\textrm{on}}$
. The inset figure shows a zoomed-in view of the 60 kpc switch-off low-powered simulations from 0 to 50 Myr showing how these small sources continue to grow as if they were active for some time. The small arrows show the time delay between the point where the jet flow stops and the point where the lobe length evolution begins to slow down.
For low-power 20 and 60 kpc switch-off simulations (see the left panel of Figure 5) there is a time delay between the jet switch off and the divergence of the length evolution from the trajectory of the active source. This is due to the travel time of the trailing jet material through the lobe. To see this more clearly, we isolate the 60 kpc switch-off simulations in the inset in Figure 5. Here, the arrows show the difference in time between the jet switch off and the point where the length evolution starts to separate from the active track. Hence, the smaller 20 and 60 kpc switch-off simulations continue to expand as if they were still active for almost 20 Myr. For the largest group simulation (Q36-v01-a25-G180; purple solid lines), we see an increase in total source expansion when the remnant phase begins. We attribute this to the low ambient density and pressure in this region of the environment (see Figure 2).
To determine when the expansion velocity of the sources matches the average sound speed in the ambient medium, we show the time derivative of the source length evolution (Figure 5) in Figure 6. To reduce the noise, we smooth the data using a LOcally WEighted Scatter plot Smoothing model (LOWESSFootnote
f
). The range of sound speeds in the undisturbed cluster and group environments are shown by the grey and red shaded regions, respectively. In the cluster environment, before
$t_{\textrm{on}}$
, the expansion rates are transonic/subsonic compared to local sound speeds. In the group, our 20 and 60 kpc switch-off sources eventually transition to subsonic velocities. We note here that these sources are strongly affected by the dynamics of the environment, making estimates of their forward velocity difficult after switch off. Simulation Q36-v01-a25-G180 continues to expand at least mildly supersonically (a Mach number between 1-2) for the duration of the simulated time.
We track the evolution of the length to width ratio of the primary lobe in Figure 7. The length to width ratio is defined as the ratio of the full length of the lobe to the maximum width measured at the lobe midpoint. We consider only the lobe propagating in the
$+z$
direction, though results are similar for the
$-z$
lobe. With this definition, low length to width ratios imply wide lobes, while high length to width ratios imply narrow ones. The initial spike in the length to width ratio tracks at
$t\sim5$
Myr corresponds to the jet punching through the central region of the host environment before the lobe has inflated. This is similar to the jet-breakout phase described in Sutherland & Bicknell (Reference Sutherland and Bicknell2007) and Turner et al. (Reference Turner, Yates-Jones, Shabala, Quici and Stewart2023).
The forward (propagating in the
$+{\mathrm{z}}$
direction) advance speed of the primary lobe shown against lobe length, illustrating the deceleration of all sources after switch-off. The figure features are consistent with Figure 5 where the left-hand panel shows all models with low injected jet kinetic powers and high powers are on the right and the thick to thin lines denote the active and remnant phases, respectively. The cross marks indicate 2
$\times t_{\textrm{on}}$
.
The time evolution of the ratio of the primary lobe length to width. The lobe width is measured as the maximum
$\pm x$
extent in a narrow slice centred at the lobe midpoint.
The switch-off time of the jets has some effect on the length to width ratio over the entire lifetime of the low-power simulations. In large, 180 kpc switch-off models, the length to width ratio stabilises, suggesting quasi-self-similar expansion after switch-off. Conversely, in 20 and 60 kpc switch-off simulations, the length to width ratio increases as trailing jet material burrows through the top of the lobe and rises away from lower radii plasma. The environment also plays a role; the higher density and pressure in the cluster environment results in narrower lobes with consistently higher length to width ratios. This occurs because the surrounding pressure in the cluster environment suppresses lateral expansion. The forward (jet-driven) expansion is less affected, allowing the jets to drill through the ambient medium.
We lastly consider the pressure profiles of all sources to probe the effect the active and remnant evolution have on the initial gas distribution. We show pressure profiles along the axis of jet propagation for all sources in Figure 8. In the active phase, the pressure profiles of low-powered sources are characteristically non-uniform, with regions of underpressure and slight overpressure. The regions of underpressure are located in the wake of the expanding head of the lobe and are more than an order of magnitude less than the pressure of the initial ambient profile. In the remnant phase, all low-powered simulations quickly return to pressure balance with the initial profile.
Pressure profiles along the axis of jet propagation for all simulations showing how remnants of low-power (top two rows) and high-power (bottom three rows) progenitors affect the initial gas distribution (solid black line). We show the pressure distributions at
$t_{on}$
(solid colored line) and then 10 Myr (dashed), 50 Myr (dot-dashed), and the maximum simulated time (indicated in Table 2, dotted line) in the remnant phase.
3.2.2. High-power progenitors
We now consider the evolution of our high-power (
$Q_{\mathrm{j}} = 10^{38}$
W), fast (
$v_{\mathrm{j}} = 0.98$
c) sources using the same methods as for the low-power sources. During the active phase, all high-power fast sources are significantly overpressured and expand supersonically through the ambient medium. The right-hand panels of Figures 5 and 6 show that simulations expanding into the group have faster rates of expansion and take roughly half the time to reach the same length; in contrast to low-power sources, where group simulations were only slightly faster and the active times for equivalent jet simulations were similar. As the lobes inflate, the length to width ratio of all high-power models decreases (see Figure 7) but will stabilise in most cases within a few Myr after switch-off, indicating quasi-self similar expansion. At later times in the remnant phase, the remnant lobe structure can become susceptible to instability growth and the weak motions of the ambient medium which cause fluctuations in the length to width ratio tracks. Since our measurement of length to width ratio takes the total lobe width at the lobe midpoint, spikes in the length to width ratio tracks of Q38-v98-a25-C180 (green lines) are reflective of ‘pinching’ that starts to develop in the main body of the lobe during the remnant phase. The onset of this ‘pinching’ can be seen in the bottom right panel of Figure 4.
The deceleration of the jet plasma and subsequent transition to a state of pressure equilibrium is set by the relative balance between jet-dominated and lobe-dominated expansion. Both wide and narrow opening angle simulations are jet-dominated within the first few Myr (as indicated by higher length to width ratios during this time). For narrow opening angles, the expansion continues to be more jet-dominated out to larger source sizes, with mildly relativistic jet velocities maintained out to approximately 90% of the lobe length. Jet-dominated simulations decelerate most rapidly once the jet flow ceases; as seen in the right-hand panel of Figure 6 however, this deceleration is also seen prior to switch-off. For wide opening angles, the transition from a jet-dominated to lobe-dominated expansion occurs quickly after the lobes begin to form.
As shown in the bottom three rows of Figure 8, all high-powered remnants are strongly overpressured in the active phase and display a uniform pressure distribution along the symmetry axis, as expected from their combination of injection power and velocity. In the remnant phase, the lobes remain overpressured (relative to the material outside the lobe/shocked ambient system) but the pressure distribution in the central region of the grid falls below the initial profile as the lobes continue to drive gas from low radii. In one simulation (Q38-v98-a7.5-C180 in the lower right of Figure 8), we see an increase in pressure of the central region at late times (the dotted green line) as gas refills the central region and pushes aside the remnant lobes.
3.3. Bow shock dynamics
We now consider the evolution of the surrounding bow shock. Initially, all simulations in this work expand supersonically, driving shocks into the ambient medium. These regions can be seen in Figures 3 and 4 as a parabolic-shaped cross-section of hot, high-density fluid surrounding each radio jet and jet-inflated lobe. We track the leading edge of the shocked region by identifying the point where the pressure gradient is maximised outside the lobe material. This matches well with the visible trace of the shocked region seen in Figures 3 and 4. We track the bow shock until the maximum simulation time is reached, the pressure gradient is less than
$1.1$
per grid cell (a necessarily low threshold to capture weak shocks while still avoiding picking up fluctuations in the ambient pressure profile), or it has exceeded the grid boundary. Low powered sources are discussed first in Section 3.3.1 and then high-power sources in Section 3.3.2.
3.3.1. Low-power progenitors
All low-power simulations produce weak bow shocks that do not considerably change in shape or separate from the lobe during the remnant phase. In the left-hand panel of Figure 9, we show the length to width ratio of the shocked region which we measure as the ratio of the forward shock radius to the lateral radius taken along the positive x-axis. Very similar results are produced if the y-axis is considered due to the spherical symmetry of the shock. For both group and cluster simulations, the length to width ratio values remain in the range of 1.5–2.2 over the entire simulated time and the shock fronts separate only marginally from the head of the lobe as shown in the right panel of Figure 9. For low-powered sources, the shock separation is more pronounced for small remnants in the group environment, likely as a result of the low-density declining atmosphere.
Evolution of the length to width ratio (left) of the shocked region, the ratio of the maximum shock front pressure to the undisturbed ambient medium (middle), the ratio of shock radius to lobe radius along the z direction (right). As before, thick lines denote active sources and thin ones denote remnants. Tracks of low-power sources are given in the top row and high-power sources are given in the bottom row. For the low-power cluster simulations, the strength of the lateral shock becomes indistinguishable from the ambient medium at around 20 Myr and so the blue lines in the right-hand panel are truncated.
The strength of the shock front is presented in the middle panel of Figure 9. At the start of the remnant phase, the pressure of the shock front for the low-power cluster simulations is, at most, twice that of the ambient medium. As a result, we are only able to track the shock front for about 20 Myr into the remnant phase before it drops below
$1.1\times P_{\textrm{env}}$
. For the group simulations, the shock pressure ratios are higher at the end of the active phase (around
$3\times P_{\textrm{env}}$
) but take over six times longer to decay to similar values as the ambient medium.
3.3.2. High-power progenitors
Our high-power, fast simulations are characterised by strong bow shocks, and they can be distinguished from the ambient medium for the total simulated time. Before the remnant phase begins, the bow shock is between 3–200 times more overpressured than the undisturbed ambient medium (see the lower middle panel of Figure 9) and is only marginally separated from the cocoon at the head of the lobe (lower right panel of Figure 9). When the remnant phase starts, the bow shock progressively separates from the head of the lobe and tends towards sphericity as shown by the decreasing length to width ratio values in the bottom left panel of Figure 9. This result is consistent with a decrease in jet head expansion velocities. A similar observation is made by Perucho et al. (Reference Perucho, Quilis and Mart2011) and (Reference Perucho, Mart, Quilis and Ricciardelli2014) for their suite of strongly relativistic jets with kinetic powers in the range of
$10^{37}$
–
$10^{39}$
W.
The separation between lobe and bow shock is most drastic in the smallest switch-off simulations, where the edge of the shocked region can increase to twice the lobe radius within 20 Myr. Meanwhile, the bow shock pressure ratio falls to below 20. The initial drop in pressure at the start of the remnant phase is most rapid for narrow opening angle simulations (the solid green track in Figure 9). For high-power, wide opening angle simulations in the group environment, the rate of bow shock pressure decrease is similar the equivalent cluster model, but it remains more than 10 times overpressured relative to the surrounding environment.
4. Discussion
In this section, we compare our results to existing numerical and analytical modelling of remnant radio sources in Section 4.1. In Section 4.2 we discuss the relative role of jet momentum flux and buoyancy in driving the expansion of our sources in the active and remnant phases and show that simple analytic expectations can provide good estimations of the onset of different evolutionary phases. In Section 4.3, we consider whether remnants can continue to provide feedback on the ambient medium after switch-off.
4.1. Benchmarking
4.1.1. Numerical modelling
We start by comparing our results to those found in similar, past works. The evolution of large-scale (several hundred kpc), high-powered (
$\gt10^{37}$
W) sources has been well studied by authors including Walg et al. (Reference Walg, Achterberg, Markoff, Keppens and Porth2014), Perucho et al. (Reference Perucho, Mart, Quilis and Ricciardelli2014) and English et al. (Reference English, Hardcastle and Krause2019). The evolution of the high-powered sources in the present study produces similar key results and features. Similar to Perucho et al. (Reference Perucho, Mart, Quilis and Ricciardelli2014) and English et al. (Reference English, Hardcastle and Krause2019) we show that the shocks and lobes tend towards sphericity once the driving forward momentum is removed. As the lobe material decelerates post switch-off, the shock front separates and weakens as it propagates into the ambient environment (see Figure 9). As expected and in line with the results of Walg et al. (Reference Walg, Achterberg, Markoff, Keppens and Porth2014), the trailing jet material disappears on timescales set by the length of the cocoon and the bulk velocity of the jet material. This is also the duration of time that the source will continue to ‘masquerade’ as if it were still active. For our high-powered sources, the jet material traverses the lobe in a shorter time than typical time between simulation outputs. Hence, these sources appear to deviate from their active behaviour immediately following switch-off.
We find the environment to be important for the dynamics of the simulated sources regardless of the jet kinetic power and injection speed. A key difference is the amount of overpressure in the lobes. For large-scale high-power sources in group environments, the source can remain significantly overpressured in the outer parts of the lobe for much longer than the previous active phase, as shown in Figure 8. This is largely attributed to the lower ambient pressure in this environment. The same behaviour is seen for low-power sources. Although the overpressure factor is lower for low-powered sources, the volume occupied by the lobes is a factor
$2-3$
larger in the group than in the cluster at a given source length. In English et al. (Reference English, Hardcastle and Krause2019), the authors find their remnant lobe dynamics to be less dependent on the environment, with all models continuing to expand at similar rates after switch-off. We show that while equivalent high-power group and cluster models decelerate at a similar rate initially, the difference in environment properties results in the lobe expansion speed of group simulations reaching the ambient sound speed later than their cluster counterparts. We further explore the transition to buoyancy in Section 4.2.
4.1.2. Analytic modelling
We now look at how our simulation results compare with analytic modelling and investigate whether existing analytic models can be used to capture the evolutionary histories of simulated remnant radio galaxies. The transition of a powerful, relativistic source from its active to remnant phase has been explicitly considered in the literature with different approximations. Krause (Reference Krause2005) derived solutions for the shock radius in a spherical thin-shell model, including the gravity of the dark matter halo. A self-similar model without gravity, but for both lobe and shock radius, was derived by Kaiser & Cotter (Reference Kaiser and Cotter2002). In that work, the authors assume a spherical radio source with uniform internal pressure and derive a set of steady-state similarity solutions that describe the evolution of the cocoon and shock radii in the active and remnant phases. In the active phase, the radio source evolution can be described by their Eq. 5. Once the jet switches off, Kaiser & Cotter (Reference Kaiser and Cotter2002) assume that the cocoon remains overpressured and immediately enters a coasting phase. For the case where the adiabatic index of the cocoon and shocked material is
$5/3$
representing non-relativistic material (a reasonable assumption for the material in the lobes of our simulations), the growth of the cocoon is given by
$R \propto t^{\alpha}$
where, for the cocoon,
$\alpha_{\textrm{c}} = \frac{2(\Gamma + 1)} {\Gamma(7+3\Gamma - 2\unicode{x03B2})} = \frac{16}{5(12-2\unicode{x03B2})}$
, and for the shocked shell,
$\alpha_{\textrm{s}} = \frac{4}{7+3\Gamma-2\unicode{x03B2}} = \frac{4}{12-2\unicode{x03B2}}$
where
$\unicode{x03B2}$
is the log slope of the ambient density profile
$\unicode{x03C1} \propto r^{-{\unicode{x03B2}}}$
.
The steady-state solutions in Kaiser & Cotter (Reference Kaiser and Cotter2002) ignore the lobe energetics during the transition between the active and remnant states. If the system takes time to settle into a coasting phase, then it will not be well described by the Kaiser & Cotter (Reference Kaiser and Cotter2002) solution. We follow the work presented in Turner et al. (Reference Turner, Yates-Jones, Shabala, Quici and Stewart2023) and Turner & Shabala (Reference Turner and Shabala2023), and derive a set of coupled differential equations that describe the expansion of the cocoon and shocked region. These are obtained from Eqs. 1, 2 and 3 in Kaiser & Cotter (Reference Kaiser and Cotter2002) but are recast such that the current expansion rate is influenced by the source’s expansion history. We assume that the radio source is elongated along the jet axis with an axis ratio, A. In this case, A is defined as the length of the lobe divided by its maximum transverse radius (not full lobe width), as given in Turner et al. (Reference Turner, Yates-Jones, Shabala, Quici and Stewart2023). The values for A that we derive from our simulations are then approximately twice those in Figure 7 at the switch-off point. Following Turner et al. (Reference Turner, Yates-Jones, Shabala, Quici and Stewart2023), we implement a standard fourth-order Runge-Kutta method which solves the following equations for cocoon and shock radii:
\begin{equation}\begin{split} \dot{R_{\textrm{c}}} & = \frac{(\Gamma - 1)Q_{\mathrm{j}}A^2 R_{\textrm{s}}^{\unicode{x03B2}}}{2\pi\Gamma R_{\textrm{c}}^2k\dot{R_{\textrm{s}}}^2} + \frac{\unicode{x03B2} R_{\textrm{c}} \dot{R_{\textrm{s}}}}{3R_{\textrm{s}}\Gamma} - \frac{2R_{\textrm{c}} \ddot{R_{\textrm{s}}}}{3\Gamma\dot{R_{\textrm{s}}}} \\[6pt] \ddot{R_{\textrm{s}}} & = \frac{3(\Gamma-1)Q_{\mathrm{j}} R_{\textrm{s}}^{\unicode{x03B2}-3}A^2}{4\pi k \dot{R_{\textrm{s}}}} + \frac{\dot{R_{\textrm{s}}}^2(2\unicode{x03B2}-3\Gamma-3)}{4R_{\textrm{s}}},\end{split} \end{equation}
where
$k = \unicode{x03C1}_{\textrm{0}} a_{\textrm{0}}^{\unicode{x03B2}}$
in which
$a_0$
and
$\unicode{x03B2}$
are, respectively, the core radius and exponent in a simple power-law description of an ambient gas density profile,
$\unicode{x03C1}(r) = \unicode{x03C1}_0(r/a_0)^{-\unicode{x03B2}}$
. Since our cosmological environments are a poor match for the simpler power-law prescription used by Kaiser & Cotter (Reference Kaiser and Cotter2002), we use the approach of Turner & Shabala (Reference Turner and Shabala2015) and approximate the simulated gas density profile as being made up of several sequential, contiguous power-laws. The
$\unicode{x03B2}$
for each is found by taking the density slope in the region of interest.
We show a comparison between a subset of our numerical simulations, the expectations of the Turner et al. (Reference Turner, Yates-Jones, Shabala, Quici and Stewart2023) approach, and the expectations of the Kaiser & Cotter (Reference Kaiser and Cotter2002) analytic equations in Figure 10. The key result is that, for sources whose active expansion is not dominated by the jet momentum thrust (e.g. our high-power, wide opening angle models), our analytical approach (orange lines) shows close agreement with the numerical simulations (grey lines). However, where the expansion is jet-dominated, both analytic approaches overestimate the source size into the remnant phase because they do not take the initial jet expansion phase into account.
A comparison of the results of the shock (dashed lines) and cocoon (solid lines) evolution for three high-power simulations. The orange and blue sets of lines respectively show the expected shock and cocoon evolution for our differential approach and the analytic solution of Kaiser & Cotter (Reference Kaiser and Cotter2002). The point at which the simulation becomes a remnant is indicated by the vertical grey dotted line.
4.2. Buoyancy
Observations from Brienza et al. (Reference Brienza2017), Mahatma et al. (Reference Mahatma2018) and Jurlin et al. (Reference Jurlin2021) show that candidate remnants can have amorphous morphologies at MHz frequencies. In fact, this type of morphology is often taken as an indicator of remnant status (Brienza et al. Reference Brienza2016). Dynamically, these relaxed structures would be expected to develop following the onset of buoyancy (e.g. Churazov et al. Reference Churazov, Bruggen, Kaiser, Bohringer and Forman2001). In this section, our goal is to determine when this takes place from our numerical results, and provide an analytical description of the buoyant regime in both the active and remnant phase. We analyse the active phase in Section 4.2.1 and the remnant phase in 4.2.2.
4.2.1. Active sources
We first must determine whether our simulations are buoyancy dominated before the onset of the remnant phase. Komissarov & Falle (Reference Komissarov and Falle1998) and Krause et al. (Reference Krause, Alexander, Riley and Hopton2012) describe a length scale,
$L_{\textrm{2}}$
, where a radio source of kinetic power
$Q_{\mathrm{j}}$
would be expected to come into pressure equilibrium with the surrounding environment,
\begin{equation} L_2 = \left(\frac{Q_{\mathrm{j}}}{\unicode{x03C1}_{\textrm{ext}} c_{\textrm{s, ext}}^3}\right)^{1/2} \end{equation}
where the external sound speed,
$c_{\textrm{s, ext}} = \left(\frac{\Gamma_{\textrm{ext}} P_{\textrm{ext}}}{\unicode{x03C1}_{\textrm{ext}}}\right)^{1/2}$
, is calculated from the local ambient pressures
$P_{\textrm{ext}}$
, and densities,
$\unicode{x03C1}_{\textrm{ext}}$
, through which the radio source propagates. We reasonably assume that the ambient medium is well-described by an ideal adiabatic index of
$\Gamma$
= 5/3. For sources with linear sizes
$L \lt L_2$
, the lobe plasma is expected to be supersonic, overpressured and the forward expansion is driven by ram pressure.
The local ambient sound speeds range from around 700–780 km/s in the cluster and 360–480 km/s in the group. Fixing the local sound speed to a median value yields
$L_2 = 16$
kpc in the cluster and
$L_2 = 70$
kpc in the group for low-power sources. The same exercise repeated for high-power sources gives
$L_2 = 160$
kpc in the cluster and 700 kpc in the group (well beyond the confines of the grid). For all low-power sources in the cluster environment, the switch-off length is in excess of the value of
$L_2$
, hence we expect these sources to be strongly buoyancy-driven prior to the start of the remnant phase. We now consider how the expected terminal rise speeds of a buoyant bubble compare to the local sound speed.
The terminal rise velocity of a buoyant bubble can be estimated following the work of Churazov et al. (Reference Churazov, Bruggen, Kaiser, Bohringer and Forman2001), Enßlin & Heinz (Reference Enßlin and Heinz2002), and Perucho et al. (Reference Perucho, Mart, Quilis and Ricciardelli2014) under the assumption that the density of the bubble is sufficiently below the density of the ambient environment. We equate the buoyant force experienced by a bubble of cross-sectional radius,
$R_{\textrm{b}}$
in an environment defined by gravity g and density
$\unicode{x03C1}_{\textrm{ext}}$
, as
$F_{\textrm{buoy}} = \frac{4}{3}\pi R_{\textrm{b}}^3g\unicode{x03C1}_{\textrm{ext}}$
, with the drag force provided by the ram pressure as the plasma moves through the external environment
$F_{\textrm{drag}} \sim C_{\textrm{d}}v_{\textrm{b}}^2\unicode{x03C1}_{\textrm{ext}} \pi R_{\textrm{b}}^2$
. The cosmological environments used in this work were selected to be as close as possible to hydrostatic equilibrium (see Paper I). In such environments, the gravitational effect,
$\unicode{x03C1}_{\textrm{ext}}(r)g$
, can be approximated as the pressure gradient in the considered region
$\frac{dP_{\textrm{ext}}}{dr} \simeq -\unicode{x03C1}_{\textrm{ext}}(r)g$
. The bubble velocity is then given by
\begin{equation} v_{\textrm{buoy}}(r)^2 = \frac{4R_{\textrm{b}} }{3C_{\textrm{d}}\unicode{x03C1}_{\textrm{ext}}(r)} \frac{dP_{\textrm{ext}}}{dr} \end{equation}
Similar to Churazov et al. (Reference Churazov, Bruggen, Kaiser, Bohringer and Forman2001) we assume a drag coefficient,
$C_{{d}} \approx 0.75$
. For a bubble of radius 3–10 kpc (the radius of our smallest sources in the active phase) moving through the inner 50 kpc of our cluster environment, the terminal rise speeds are expected to be half of the local sound speed. We expect that these sources are still significantly driven by pressure during the active phase. For larger sources, the terminal rise velocity in the region 50–90 kpc from the environment core is comparable to the sound speed. Based on our estimate of
$L_2$
and the expected rise velocities of sources in this region, we do not expect any high-power simulations to be buoyant before the remnant phase.
4.2.2. Remnant sources
Once the jet injection into the cocoon stops and
$Q_{\mathrm{j}}$
is set to 0, a new length scale needs to be considered. Using the same approach as Komissarov & Falle (Reference Komissarov and Falle1998) and Krause et al. (Reference Krause, Alexander, Riley and Hopton2012), we derive the expected lobe length for a remnant radio source terminated at length
$L_{\textrm{on}}$
and time
$t_{\textrm{on}}$
by considering where the lobe expansion speed is equal to the external sound speed,
\begin{equation}L_{2\textrm{, rem}} = L_{\textrm{on}} \left(\frac{c_{\textrm{s, ext}} t_{\textrm{on}} }{\alpha L_{\textrm{on}}} \right)^{\frac{\alpha}{\alpha-1}} \end{equation}
where
$\alpha$
is the slope of the temporal length evolution for a particular source,
$L \propto t^{\alpha}$
, and is taken directly from the simulations. Generally, we find a very good agreement between the majority of our simulations and the calculated
$L_{2\textrm{, rem}}$
value, as shown in Figure 11. For high-power simulations in the cluster, the calculated
$L_{2\textrm{, rem}}$
value agrees remarkably well with the numerical results (see the orange and green lines in the top panel of Figure 11). For those simulations that do not reach the sound speed during the simulated time (i.e. our large switch-off simulations in the group environment; see the red and purple solid lines in the bottom panel of Figure 11), the
$L_{2, \textrm{rem}}$
value is larger than the computational domain. For these simulations, we predict
$L_{2, \textrm{rem}}$
values of 350 and 790 kpc for Q38-v98-a25-G180 and Q36-v01-a25-G180 respectively.
Forward advance velocities for a subset of our simulations showing the active (thick) and remnant (thin line) phases for these simulations. The predictions of the buoyant, remnant length scale
$L_{2,\textrm{rem}}$
are denoted by the triangle markers. The top panel shows cluster simulations and the bottom panel shows group simulations. For the large group simulations shown in the lower panel, the remnant length scale is larger than the grid size and hence, triangle markers are not shown. The cross marks indicate
$2t_{\textrm{on}}$
and the grey shaded regions show the range of sound speeds in the ambient medium.
The median tracer value as a function of lobe length along the axis of jet propagation (such that material at the lobe head has a lobe length value of 1) for our largest low-powered simulation in the cluster (Q36-v01-a25-C180, left) and largest high-powered simulation in the cluster (right). We show the tracer distribution at four time steps: at switch-off (yellow), and at 5, 40, and 80 Myr into the remnant phase (pink, purple, black, respectively). The slow changes in tracer values for low-powered sources suggests mixing occurs slowly, while tracer values in the high-powered simulation drop rapidly, indicating fast mixing.
Using Eq. 8 we again consider the terminal rise speeds of our simulations during the remnant phase. For our simulations in the cluster environment, the forward velocities for 20 and 60 kpc switch-off sources fall below the 300–350 km/s buoyant velocity in the inner 50 kpc of the cluster profile during the remnant phase. The result is similar for 20 and 60 kpc simulations in the group environment where rise speeds are
$\lesssim 250$
km/s. This is consistent with the estimated terminal rise speeds of small bubbles outlined in Section 4.2.1. Remnant lobes are expected to rise at a velocity below the buoyant velocity once their density approaches that of the surrounding environment, continuing to decelerate as drag becomes significant (Churazov et al. Reference Churazov, Forman, Jones and Bohringer2000). This is true of our low-powered sources, where the density of the lobe material becomes increasingly similar to the nearby ambient medium as the lobes mix and entrain surrounding material (we refer the reader to the top two rows of Figure 4). For high-powered sources, the lobe density is instead a few orders of magnitude lower than the ambient medium at switch-off, however, as these sources remain as a single mass (as opposed to two detached lobes which could independently rise apart) centred at the gravitational potential, buoyancy has a greatly reduced effect on the lobe forward velocities. Instead, the sources drift off the axis of jet propagation.
For larger sources which reach 100 kpc or more, Rayleigh-Taylor instabilities begin to pinch the lobe material, enabling dense ambient medium to fall inward towards the gravitational potential. Buoyancy then becomes more important in driving the lobes away from the centre of the grid.
The results obtained in this work can be compared to the similar analysis conducted in Perucho et al. (Reference Perucho, Mart, Quilis and Ricciardelli2014) for their suite of two-dimensional axisymmetric simulations spanning jet powers in the range of
$10^{37}$
–
$10^{39}$
W. For their high-powered simulations (e.g. their model J46; a
$10^{39}$
W,
$v_{\mathrm{j}}=0.9$
c jet), they find the measured advance speed of the cavity to exceed the predicted terminal rise speed by a factor of two; too high to be explained by buoyant rise alone. Conversely for their low-powered sources the measured rise speed sits below the predicted value indicating buoyancy is more likely. In this work, we find strong evidence for the buoyancy-driven expansion of low-powered sources even before the jets have switched off in some cases (our large, low-powered sources). However, for our high-powered sources, we do find terminal rise speeds that match the measured forward advance speeds of the lobe in the case where the environment is more dense (our cluster simulations).
4.3. Heating from remnant lobes and shocks?
In this section, we discuss the ability of remnant lobes and shocks to continue to heat the surrounding medium in the remnant phase. This is particularly important as the lobes and shocked regions of our remnant sources tend to sphericity after the jets switch off, offering a mechanism to heat the surrounding medium isotropically. This is feature not seen in active sources due to the inherent directionality of their expansion; see also Vernaleo & Reynolds (Reference Vernaleo and Reynolds2006).
The ratio of bow shock front temperature to environment temperature for all simulations during the active (thick lines) to the remnant phase (thin lines).
Feedback from AGN can act as a thermostat for the host environment, quenching otherwise catastrophic cooling at the centres of clusters (Yang & Reynolds Reference Yang and Reynolds2016). The channels through which this heating may occur include cavity heating (Churazov et al. Reference Churazov, Bruggen, Kaiser, Bohringer and Forman2001), the driving of strong shocks (Fabian et al. Reference Fabian2003; Perucho et al. Reference Perucho, Mart, Quilis and Ricciardelli2014; Shabala, Kaviraj, & Silk Reference Shabala, Kaviraj and Silk2011), the dissipation of sound waves, the mixing of thermal gas within the jet-inflated remnant lobes (Hillel & Soker Reference Hillel and Soker2016), and the uplifting of low-entropy gas to higher radii by buoyantly rising remnant lobes (Chen et al. Reference Chen, Heinz and Enßlin2019). Shabala et al. (Reference Shabala, Kaviraj and Silk2011) and Perucho et al. (Reference Perucho, Mart, Quilis and Ricciardelli2014) have shown that the feedback mechanism is dependent on the initial power of the jets. High-powered sources are responsible for the driving of strong shocks, while low-powered sources heat the environment more gently by way of weak shocks and the mixing of the lobes with the ambient medium. Regarding the driver of this mixing, Bourne & Sijacki (Reference Bourne and Sijacki2021) have suggested that cluster weather, rather than the development of Kelvin-Helmholtz and Rayleigh-Taylor instabilities, drives lobe evolution after the lobe-inflation phase. As low-powered jet lobes are more readily disrupted by the cluster weather, they more readily mix with the ICM. On the other hand, Hillel & Soker (Reference Hillel and Soker2016, Reference Hillel and Soker2017) have argued that it is the development of instabilities that drive mixing between the shocked ambient medium and the lobes.
The passive tracer advected with the jet fluid is a good indicator of the composition of jet material and ambient material in the lobes of our simulated sources. Considering the change in tracer value as a function of length over time gives an indication of where the mixing is taking place and how rapidly it occurs after the jet switches off. At injection, jet fluid is assigned a tracer value of 1 while the initial ambient medium is assigned 0. Tracer values in between 0 and 1 indicate mixed ambient and jet material. In Figure 12 we consider the median tracer value within each row of grid cells along the axis of jet propagation (
$z$
in our set-up) for two representative simulations; our largest low-powered cluster simulation (left panel of Figure 12) and our largest high-powered cluster simulation (right panel). To isolate the lobe material, we have excluded the recently injected jet material by removing all material with a tracer value greater than 0.2 and above the tracer cutoff of
$10^{-4}$
.
For low-powered sources, the tracer distributions evolve slowly over time and mixing appears to occur preferentially at the base of the lobe. For high-powered sources, tracer values drop rapidly after switch-off, particularly at the head of the lobe. The contact surface is Rayleigh-Taylor unstable due to the deceleration of the denser shocked gas shell with respect to the cocoon (Gull & Northover Reference Gull and Northover1973). The instability growth timescale is proportional to the density contrast, which is 2–3 orders of magnitude higher for the more powerful jet. A more detailed study of the post-switch off feedback channels for remnant lobes of high- and low-powered progenitors is deferred to later work.
We next consider the heating by the bow shock on the ambient medium. In Figure 13, we present the temporal evolution of the ratio of shock temperature to the environment temperature for low- (left panel) and high-powered (right panel) simulations. As with Figure 9 in Section 3.3, we take the measurement of shock temperature as the maximum value at the leading edge of the shocked region along the propagation axis. For some low-powered simulations (namely, the largest low-power simulations Q36-v01-a25-C180 and Q36-v01-a25-G180), the temperature ratio rises after switch-off, indicating the shock temperature decreases more slowly than the ambient profile. The temperature ratio for high-powered simulations is consistently higher than low-powered simulations, supporting the idea that the driving of strong shocks is an important heating mechanism for high-powered sources. Further to this, the temperature ratio of high-powered objects plateaus after an initial sharp drop at the start of the remnant phase. Given the significant separation of the shock from the lobe material presented in Section 3.3, the shocks of high-powered sources present a way to isotropically heat the surrounding environment at distances up to twice the lobe radii well into the remnant phase. How efficient the feedback is, depends on the evolution of the post-shock gas. For example, relatively flat environments and non-negligible thermal conductivity can both lead to slow cooling of the post-shock gas and hence efficient AGN feedback (Binney & Tabor Reference Binney and Tabor1995; Alexander Reference Alexander2002). Overall, our results support the idea that both feedback channels are viable. In Raouf et al. (Reference Raouf, Shabala, Croton, Khosroshahi and Bernyk2017), both ‘shock’ and ‘bubble’ modes of RLAGN feedback are considered in a semi-analytic galaxy formation model. They show both are needed to explain the observed galaxy properties, namely the jet-power luminosity relation, stellar mass-radio luminosity relation, and the radio luminosity function as it scales with redshift.
5. Conclusions
We have performed three-dimensional numerical hydrodynamic simulations of relativistic and sub-relativistic AGN jets propagating into realistic environments extracted from cosmological simulations. We evolve these sources from the active phase into the remnant phase by terminating the flow of injected material onto the simulation grid. As an extension to previous works in this area, we consider a parameter space of both low- (Q
$_{\mathrm{j}} = 10^{36}$
W) and high-power (Q
$_{\mathrm{j}} = 10^{38}$
W) jets of varying sizes, in two different environments. The broad dynamics and behaviour of our simulated sources are consistent with earlier work (e.g. Walg et al. Reference Walg, Achterberg, Markoff, Keppens and Porth2014; Perucho et al. Reference Perucho, Mart, Quilis and Ricciardelli2014; English et al. Reference English, Hardcastle and Krause2019) and can be understood analytically by combining the approaches of Kaiser & Cotter (Reference Kaiser and Cotter2002) and Turner & Shabala (Reference Turner and Shabala2015). We find the following key results:
-
(i) We confirm several prior results about the dynamic evolution of remnant radio galaxies. We find that some radio remnants from high-power progenitors tend towards more spherical shapes while the leading shock progressively separates from the lobe material. We also confirm the expectation of Walg et al. (Reference Walg, Achterberg, Markoff, Keppens and Porth2014) that remnants will experience a delay in their dynamic transition from active to remnant sources. This delay period is set by the time it takes for the trailing jet material to traverse the lobe. For low-powered, slow simulations, this can be as long as 5–10 Myr.
-
(ii) We validate the use of analytic modelling approaches to track the evolution of some remnant radio sources. The lobe and shock evolution of our high-powered, wide jet half-opening angle remnants is well described by analytic approaches if the expansion history is taken into account, allowing for the lobe deceleration and shock separation to be captured correctly. For sources with a narrow jet half-opening angle, the expansion is primarily driven by the momentum of the jet. After switch off, lobe expansion decelerates rapidly, which is not captured in analytic models such as Kaiser & Cotter (Reference Kaiser and Cotter2002).
-
(iii) We show that analytic approaches can be used effectively to capture the transition from the pressure driven, to the buoyant phase of some simulated remnants. We predict that the overpressure in the lobes of remnants is sufficient to extend the pressure-driven phase by over 100 Myr post switch-off in low density environments.
-
(iv) Remnants of high-power progenitors can supply significant feedback to the ambient gas through bow shock heating, which continues into the remnant phase. We speculate that the feedback supplied by low-powered sources is inefficient, with mixing between the shock and lobe occurring very gradually.
Detailed predictions of observable synchrotron emission are required to interpret observations of remnant candidates, and complement the dynamical study presented here. We defer this to the companion work in Stewart et al. (in preparation).
Acknowledgements
We thank the anonymous referee for their helpful comments which added value to the manuscript. GS thanks the Australian Government for an Australian Government Research Training Program RTP scholarship and the CSIRO for an ATNF studentship. SS, CP and PYJ acknowledge the Australian Research Council grant DP240102970. This research was carried out using high-performance computing infrastructure provided by the Tasmanian Partnership for Advanced Computing (TPAC). We acknowledge Geli Kourakis and John Miezitis who provided valuable assistance throughout this project. We acknowledge the support of the developers providing the Python packages that made this work possible: JupyterLab (Kluyver et al. Reference Kluyver2016), Matplotlib (Hunter Reference Hunter2007), NumPy (Harris et al. Reference Harris2020), SciPy (Virtanen et al. Reference Virtanen2020), and Astropy (Robitaille et al. Reference Robitaille2013).
Open access
