# Numerical modelling of the lobes of radio galaxies - Paper V: Universal
Pressure Profile cluster atmospheres

M. Stimpson,1 M.J. Hardcastle,1 M. G. H. Krause1
1Centre for Astrophysics Research, Department of Physics, Astronomy and
Mathematics, University of Hertfordshire, College Lane, Hatfield,
Hertfordshire AL10 9AB, UK
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present relativistic magnetohydrodynamic modelling of jets running into
hydrostatic, spherically symmetric cluster atmospheres. For the first time in
a numerical simulation, we present model cluster atmospheres based upon the
Universal Pressure Profile (UPP), incorporating a temperature profile for a
‘typical’ self-similar atmosphere described by only one parameter - $M_{500}$.
We explore a comprehensive range of realistic atmospheres and jet powers and
derive dynamic, energetic and polarimetric data which provide insight into
what we should expect of future high-resolution studies of AGN outflows. From
the simulated synchrotron emission maps which include Doppler beaming we find
sidedness distributions that agree well with observations. We replicated a
number of findings from our previous work, such as higher power jets inflating
larger aspect-ratio lobes and the cluster environment impacting the
distribution of energy between the lobe and shocked regions. Comparing UPP and
$\beta$-profiles we find that the cluster model chosen results in a different
morphology for the resultant lobes with the UPP more able to clear lobe
material from the core; and that these different atmospheres influence the
ratio between the various forms of energy in the fully developed lobes. This
work also highlights the key role played by Kelvin-Helmholtz (KH)
instabilities in the formation of realistic lobe aspect-ratios. Our
simulations point to the need for additional lobe-widening mechanisms at high
jet powers, for example jet precession. Given that the UPP is our most
representative general cluster atmosphere, these numerical simulations
represent the most realistic models yet for spherically symmetric atmospheres.

###### keywords:

galaxies: jets, clusters – methods: numerical

††pubyear: 2022††pagerange: Numerical modelling of the lobes of radio galaxies
- Paper V: Universal Pressure Profile cluster atmospheres–B

##  1 Introduction

The modelling of the lobes of AGN jets began with analytical approaches such
as the key work of Scheuer (1974) and later on Kaiser & Alexander (1997).
Access to high-performance computing in recent years has fuelled modelling
using numerical simulations: varying the input parameters of our models and
comparing the output with observations tells us much about the conditions and
processes in the AGN. Early simulations, such as the 2DHD model of Norman et
al. (1982) confirmed the predictions made by Blandford & Rees (1974) about
bipolar jets creating hot spots and a bow shock at the jet head, followed by a
contact discontinuity between the lobe and shocked ambient medium. More recent
models incorporate more of the physics of actual AGN’s, such as the 3DRMHD
work of Mignone et al. (2010) who find that 3D models reveal kink
instabilities which are not seen in similar 2D models. Recent reviews of
numerical modelling of jets include Martí (2019), Komissarov & Porth (2021)
and Bourne & Yang (2023). A key role for modelling now is to predict the
images which will be captured by the next generation of X-ray satellites such
as Athena (Nandra et al., 2013), and of radio telescopes such as the Square
Kilometer Array (Carilli & Rawlings, 2004).

The primary factor which shapes large-scale jet structure is the power of the
jet (Rawlings & Saunders, 1991; Worrall, 2009; O’Dea & Saikia, 2021). Using a
method based upon observations of the jet terminal hotspots, Godfrey & Shabala
(2013) obtained jet powers for FR IIs in the range ${10}^{38}$ to
${10}^{39}$W; in good agreement with other researchers using a variety of
techniques (Ghisellini et al., 2014; Gitti et al., 2012). Kaiser & Alexander
(1997) found that greater powers ($Q>${10}^{37}$$W) formed FR IIs; similar
relations have been demonstrated by the simulations of Massaglia et al. (2016)
and Ehlert et al. (2018). Furthermore, relativistic MHD simulations conclude
that at the intermediate power between FR I and FR II, it is other factors
such as the density ratio, Lorentz factor and magnetic field strength which
determines whether the lobe is FR I or FR II (Mignone et al., 2010; Mukherjee
et al., 2020; Massaglia et al., 2022). Massaglia et al. (2022) note that the
observation epoch also plays a role as many features are time-dependent.
Perpendicular shocks close to the radio core have been shown to form in jets
with wide initial opening angle (Krause et al., 2012; Yates-Jones et al.,
2022a). Particle acceleration at such shocks could be responsible for the
brightening of FR I jets at flaring points.

The structure of the cluster atmosphere is also a key factor which will
influence the morphology of the lobes. Early researchers used uniform
distributions (Norman et al., 1982; Koessl & Mueller, 1988; Lind et al., 1989)
which could only represent real atmospheres over short distances; an
improvement on these is the now widely-used and more realistic $\beta$-model
(Reynolds et al., 2002; Basson & Alexander, 2003; Zanni et al., 2003; Krause,
2005). Observations of powerful radio sources has demonstrated that
asymmetries in the distribution of ionized gas is correlated with the
structural asymmetry of the radio lobes (Pedelty et al., 1989; McCarthy et
al., 1991; Gopal-Krishna & Wiita, 2000), which indicates that environmental
asymmetries play a role in creating structural asymmetries in the radio lobes;
as demonstrated in numerical models of jets propagating through inhomogeneous
environments (e.g. Jeyakumar et al., 2005; Gaibler et al., 2009, 2011; Yates-
Jones et al., 2021, 2022b); in particular, Tanner & Weaver (2022) demonstrated
that lower power jets are impacted more and Wagner & Bicknell (2011) and
Wagner et al. (2012) showed that inhomogeneities in the ICM impact the
transfer of energy from the lobe to the surroundings. Furthermore, using an
environment derived from a simulation of a dynamically active cluster
represents another step towards greater realism (e.g. Heinz et al., 2006;
Mendygral et al., 2012). Such large scale motion of the ICM disrupts the jets
and plays a significant role in spreading out the energy injected (Morsony et
al., 2010; Bourne & Sijacki, 2017; Bourne et al., 2019; Bourne & Sijacki,
2021). Dynamic modelling has also been used to investigate self-regulated
feedback (Meece et al., 2017; Ehlert et al., 2023), accretion rates (Prasad et
al., 2015) and energy transfer mechanisms Ehlert et al. (2018). For reviews
see Meece et al. (2017) and Bourne & Yang (2023). In this work we model a
spherically symmetric hydrostatic profile; similar models have been used by
others to investigate feedback mechanisms and energy transfer (e.g. Yang &
Reynolds, 2016b, a; Weinberger et al., 2017; Weinberger et al., 2023).
However, for the first time in a numerical simulation, we use the Universal
Pressure Profile (Arnaud et al., 2010) as a generalised hydrostatic,
spherically symmetric cluster atmosphere.

The numerical model used in this study is a development of that described in
Hardcastle & Krause (2013), Hardcastle & Krause (2014), English et al. (2016)
and English et al. (2019) (henceforth referred to as ; ;  and Paper 4). In all
these previous papers we used a $\beta$-atmosphere, but here we use the more
realistic UPP atmosphere, as well as a jet with a higher Lorentz factor
($\gamma=10$) which matches well with the values seen on parsec scales. The
model also employs stretched grids in order to model the central regions at a
much higher resolution and so enable a considerably narrower, and so more
realistic, injection cylinder. In Section 2 we will present the UPP atmosphere
used in this study, including the temperature profile and the ambient magnetic
field. In Section 3 we describe how the simulated atmosphere is implemented
and the jet parameters used. Our results are presented in Section 4; in
Section 5 we discuss the findings and our Summary and Conclusions are found in
Section 6.

##  2 Cluster Atmospheres

###  2.1 Galaxy cluster density and pressure distributions

The X-ray emission from the hot plasma in galaxy clusters was studied by
Cavaliere & Fusco-Femiano (1976) who developed a model in which the material
of the cluster has a density profile described by King’s approximation (King,
1962). This self-consistent isothermal model is referred to as the
$\beta$-model and is widely used to describe the density profile in clusters
of galaxies, it can be written

$\rho(r)=\rho_{0}\left[1+\left(\frac{r}{r_{c}}\right)^{2}\right]^{-\frac{3\beta}{2}}$
(1)

where $\rho_{0}$ is the density at the centre of the cluster, $r_{c}$ is the
‘core radius’ and $\beta$ is an indication of the gradient beyond the core
radius. Despite its wide use in simulations, it has long been recognised that
the $\beta$-model does not adhere well to observations, particularly near the
core where observations indicate a ‘cusp’ rather than the constant density
produced by the beta model (Frenk et al., 1988).

The Navarro, Frenk and White (NFW) model also describes the density profile of
the cluster and is an improvement on the $\beta$-model. The authors describe
the NFW as a ‘universal’ density profile; it was the result of conducting a
series of N-body simulations of higher resolution than previous studies
(Navarro et al., 1995). This removed the constant density of the $\beta$-model
and replaced it with more of a ‘cusp’.

Based upon observations of X-ray clusters with Chandra and on numerical
simulations on scales larger than these, Nagai et al. (2007) proposed a
‘generalised’ NFW (GNFW) model in which they parametrized the profile further;
instead of expressing their distribution in terms of density (as the $\beta$
and the NFW models were) they expressed it in terms of the gas pressure of the
cluster. The version below is that presented by Arnaud et al. (2010) and is
written in terms of the average scaled pressure $\mathbbm{p}$ at a normalized
distance $x$ (Equation 3) from the cluster centre, the profile is

$\mathbbm{p}(x)=\frac{P_{0}}{\left(c_{500}x\right)^{\gamma}\left[1+\left(c_{500}x\right)^{\alpha}\right]^{\left(\beta-\gamma\right)/\alpha}},$
(2)

where $P_{0}$ is the pressure at the centre of the cluster and the parameters
$\gamma,\alpha,\beta$ are respectively the central slope $\left(r\ll
r_{s}\right)$, intermediate slope $\left(r\sim r_{s}\right)$ and outer slope
$\left(r\gg r_{s}\right)$. The scale radius $r_{s}$, is defined as the radius
where the logarithmic slope of the density profile is $\alpha=-2$; and the
concentration is defined as $c_{500}\equiv R_{500}/r_{s}$. $R_{500}$
represents the radius of the cluster corresponding to a mean mass density
contrast of $500$ times the critical density of the Universe. These
parametrized values are linked to real values of pressure $P(r)$ and radial
distance $r$ using the scaling relations

$P(r)=p(x)P_{500}\qquad\text{and}\qquad x\equiv r/R_{500}$ (3)

where $p(x)$ is the normalized pressure and is linked to the average scaled
profile $\mathbbm{p}(x)$ by an empirical term which reflects the deviation
from standard self-similar scaling:

$p(x)=\mathbbm{p}(x)\left[\frac{M_{500}}{$3\text{\times}{10}^{14}$\text{h}^{-1}_{70}M_{\odot}}\right]^{\alpha(x)}$
(4)

where $\alpha(x)$ is a variable in $x$ linked to the mass of the cluster and
the dimensionless Hubble constant h${}_{70}=$ h/H0 where
H${}_{0}=70\,$kms-1Mpc-1. $P_{500}$ is the ‘characteristic pressure’ which is
dependent upon mass and redshift as follows

$P_{500}=$1.65\text{\times}{10}^{-3}$h(z)^{8/3}\left[\frac{M_{500}}{$3\text{\times}{10}^{14}$\text{h}^{-1}_{70}M_{\odot}}\right]^{2/3}\text{h}^{2}_{70}\,\text{keV}\,\text{cm}^{-3}$
(5)

where $h(z)$ is the ratio of the Hubble constant at redshift $z$ to its
present value, $h(z)=H(z)/H_{0}$. $M_{500}$ is the mass contained within the
radius $R_{500}$ at which the mean mass density is $500$ times that of the
critical density of the Universe at the cluster redshift $\rho_{c}(z)$.
$M_{500}$ and $R_{500}$ can be found from one another; from the definition of
$M_{500}$ we have

$M_{500}=\frac{4\pi}{3}R^{3}_{500}500\rho_{c}(z)\qquad\text{where}\qquad\rho_{c}(z)=\frac{3H(z)^{2}}{8\pi
G}$ (6)

where $G$ is the gravitational constant and
$H(z)=H_{0}\sqrt{\Omega_{M}(1+z)^{3}+\Omega_{\Lambda}}$ where, for a flat
$\Lambda$CDM cosmology, $\Omega_{M}=0.3$ and $\Omega_{\Lambda}=0.7$. These
relations are the GNFW model and the correct choice of parameters will result
in a very good fit to the pressure profiles of galaxy clusters (as shown in
Appendix C of Arnaud et al. (2010)). From the UPP cluster pressure profile
described here and using a cluster temperature profile (see Section 2.3) we
can recover the density profile using the ideal gas law $p=nkT$. A comparison
between the $\beta$ and the UPP density profiles can be seen in Figure 11.

###  2.2 The Universal Pressure Profile

Table 1: UPP parameters. The data used in Arnaud et al. (2010) is based upon 33 local ($z<0.2$) clusters drawn from the REFLEX catalogue and observed with XMM-Newton with mass in the range $${10}^{14}$M_{\odot}<${10}^{15}$M_{\odot}$. The data used by He et al. (2021) is derived from simulations compared with X-ray data (REXCESS). References | $P_{0}$ | $c_{500}$ | $\gamma$ | $\alpha$ | $\beta$
---|---|---|---|---|---
Arnaud et al. (2010) | 8.130$\text{h}_{70}^{-3/2}$ | 1.156 | 0.3292 | 1.0620 | 5.4807
He et al. (2021) | 5.048$\text{h}_{70}^{-3/2}$ | 1.217 | 0.433 | 1.192 | 5.490

Arnaud et al. (2010) derived an average GNFW profile (see Table 1 for the
parameter values). For their choice of parameters the scaled pressure profiles
do not show any significant dependence on mass; in other words Equation 4
reduces to $p(x)=\mathbbm{p}(x)$ and so their model is self-similar.

When X-ray measurements are used to estimate cluster masses, it is assumed
that the cluster is in equilibrium (i.e. a perfectly relaxed cluster) (e.g.
Arnaud et al., 2010); however, observations indicate that clusters are not all
relaxed as non-thermal pressure support is also present (e.g. Schuecker et
al., 2004; Sanders et al., 2011; Walker et al., 2015; Hitomi Collaboration et
al., 2016; Hofmann et al., 2016; Eckert et al., 2019; Siegel et al., 2018);
furthermore, simulations of clusters undergoing mergers or feedback processes
substantiate these observations (e.g Vazza et al., 2012; Vazza et al., 2017;
Nelson et al., 2014; Gupta et al., 2017; Bennett & Sijacki, 2022)
demonstrating that cluster formation leads to significant non-thermal gas
processes such as turbulent flows and bulk motions. Neglecting the kinematics
leads to a systematic underestimation of the masses of galaxy clusters: this
is the hydrostatic mass bias. He et al. (2021) describe how they employ a
simulation (the Mock-X analysis framework) devised by Barnes et al. (2021)
which is able to model the evolution of clusters, including the non-thermal
pressure support, and simulate X-ray emission. Their study leads to the
debiased values for the GNFW parameters (see Table 1) which can be considered
to be more accurate than those provided by previous studies; in addition, they
confirmed the self-similarity conclusion of Arnaud et al. (2010) in that this
set of generalised parameters does not depend upon mass. Their findings
indicate that the UPP of Arnaud et al. (2010) is $\approx 5\%$ higher than
their debiased pressure values at the centre, rises to $\approx 20\%$ at
$R_{500}$ and reaches almost $\approx 35\%$ in the outermost regions. We will
use the debiased UPP values of He et al. (2021) in this study.

###  2.3 Cluster temperature profile and magnetic fields

In a cool core cluster the temperature rises steeply away from the centre and
reaches a peak around a tenth of $R_{500}$, then reduces gradually towards
large radii. The core is believed to be at a lower temperature as a result of
radiative cooling; the inner regions are at a higher pressure than the
$\beta$-model predicts and this leads to a greater luminosity and so shorter
cooling time than the rest of the cluster. Without compensation through a
heating mechanism, the core temperature falls (Fabian, 1994). Vikhlinin et al.
(2006), using Chandra data for a sample of nearby relaxed clusters, derived
the following radial variation of temperature $T(r)$:

$\frac{T(r)}{T_{mg}}=1.35\frac{(x/0.045)^{1.9}+0.45}{(x/0.045)^{1.9}+1}\frac{1}{(1+(x/0.6)^{2})^{0.45}}$
(7)

where $x\equiv r/R_{500}$ is the parametrised radial distance and $T_{mg}$ is
the gas-mass-weighted temperature, defined by Vikhlinin et al. (2006) as:

$R_{500}hh(z)=830\left(\frac{T_{mg}}{5\text{keV}}\right)^{1.47/3}$ (8)

where $h=0.72$ and all other quantities are defined above. The authors point
out that their temperature profiles are self-similar when scaled to the same
overdensity radius, which is in agreement with previous authors. In addition,
they cite good agreement between their $M-T$ relation and that produced by
other authors (e.g Arnaud et al., 2005). In this study we will employ the
temperature profile described by Vikhlinin et al. (2006).

The theory of cluster magnetic fields suggests that they scale as $B\propto
n_{e}^{1/2}$; further details can be found in Appendix A.1.

##  3 SIMULATION SETUP

###  3.1 Creating a Universal Pressure Profile (UPP) model atmosphere

The UPP is a self-similar profile with one input variable: cluster mass111We
assume that the model cluster is at a distance of $z=0$.; it represents
observed cluster profiles more faithfully than any other model atmosphere. The
cluster mass is implemented as $M_{500}$ in units of
$\text{h}^{-1}_{70}\times${10}^{14}$M_{\odot}$. From observations Planck
Collaboration et al. (2013) provide values for $M_{500}$ for clusters in the
range $1$ to $$15\text{\times}{10}^{14}$\,\text{h}^{-1}_{70}M_{\odot}$;
although other authors have found that the UPP is a good fit outside this
range (e.g. Sun et al. (2011) in their study of groups of galaxies in the
range ${10}^{13}$ to
$$1.5\text{\times}{10}^{14}$\,\text{h}^{-1}_{70}M_{\odot}$); examples of
clusters towards the upper limit are very rare. We will use $M_{500}$ values
in the range $0.33$ to
$$9\text{\times}{10}^{14}$\,\text{h}^{-1}_{70}M_{\odot}$; the parameters are
summarised in Table 2. Pressure, temperature, density and entropy
distributions for the UPP cluster atmosphere for a range of values of
$M_{500}$ are displayed in Figure 1. These average pressure profiles can be
compared with those of individual clusters from the REXCESS sample, upon which
they are based (Arnaud et al., 2010; He et al., 2021); as well as temperature
profiles (Vikhlinin et al., 2006) and density profiles (Croston et al., 2008).
Furthermore, the variation of entropy with distance is derived from the
temperature and density values (see caption) and compares favourably with
profiles derived from observations (e.g. Donahue et al., 2006; Ghirardini et
al., 2017; Babyk et al., 2018).

Table 2: UPP and equivalent $\beta$-profile atmosphere parameters used to create the pressure, temperature and density distributions for a range of realistic cluster atmospheres. The $M_{500}$ value represents the mass of the cluster; the other six values are input into PLUTO. $P_{500}$, $R_{500}$, $p_{0}$ and $r_{c}$ are all input in simulation units of pressure, length, pressure and length; whereas $T_{mg}$ is in Kelvin. The parameters $p_{0}$, $r_{c}$ and $\beta$ define the equivalent $\beta$-profile (see Section 5). Each run has a label (e.g. jetXX_haloYY) where XX represents the power of the jet with values of $0.5,1,2$ or $4\text{\times}{10}^{38}$W; which are indicated by XX having values of $05,10,20$ and $40$ (see Table 3) and YY indicates the mass as shown in this table. run | $M_{500}$ | $P_{500}$ | $R_{500}$ | $T_{mg}$ | $p_{0}$ | $r_{c}$ | $\beta$
---|---|---|---|---|---|---|---
| ($\times${10}^{14}$\text{h}^{-1}_{70}M_{\odot}$) | (sim.) | (sim.) | (K) | (sim.) | (sim.) |
jetXX_halo03 | 0.3333 | $2.2598\text{\times}{10}^{-8}$ | 232.90 | $1.0076\text{\times}{10}^{7}$ | $5.038\,91\text{\times}{10}^{-7}$ | 16 | 0.48
jetXX_halo10 | 1 | $4.7007\text{\times}{10}^{-8}$ | 335.91 | $2.1276\text{\times}{10}^{7}$ | $9.2651\text{\times}{10}^{-7}$ | 27 | 0.51
jetXX_halo30 | 3 | $9.7778\text{\times}{10}^{-8}$ | 484.46 | $4.4922\text{\times}{10}^{7}$ | $1.9059\text{\times}{10}^{-6}$ | 39 | 0.50
jetXX_halo90 | 9 | $2.0339\text{\times}{10}^{-7}$ | 698.71 | $9.4849\text{\times}{10}^{7}$ | $3.9827\text{\times}{10}^{-6}$ | 72 | 0.70
Figure 1: UPP pressure (top), temperature, density and entropy (bottom)
profiles for the values of $M_{500}$ (in units of
$\times${10}^{14}$\text{h}^{-1}_{70}M_{\odot}$) used in this study. Distance
is from the centre of the cluster. Pressure and density are measured in
simulation units, in which $p_{0}=$2.7\text{\times}{10}^{-6}$\,$Pa and
$\rho_{0}=$3.01\text{\times}{10}^{-26}$\,$g/cm3. Specific entropy is
calculated as $K=Tn_{e}^{-2/3}$ where $T$ is the temperature and $n_{e}$ is
the electron number density (Voit et al., 2005); here we show values
normalised to the value at a distance of $R_{500}$. The gradual decrease in
temperature with distance beyond the maximum is not easily seen for this range
of distance, a more significant decrease would be visible out to $\sim\,$1
Mpc.

Following the methodology of Huarte-Espinosa et al. (2011), Paper 2, Paper 3
and Paper 4; we implement a cluster magnetic field which is multi-scaled,
tangled and has a magnitude related to the cluster density profile. A summary
of the techniques used to create this magnetic field can be found in Appendix
A.2. The pressure and density profiles are interpolated, along with the three
magnetic vector potentials, into the PLUTO domain at the initialisation step.

###  3.2 Ensuring the gravitational stability of the UPP model atmosphere

Starting with the equation for hydrostatic equilibrium for a spherically
symmetric cluster, the pressure exerted by the atmosphere $p(r)$ is assumed to
be balanced by the dark matter potential $\Phi(r)$ as

$\frac{dp(r)}{dr}=-\rho(r)\frac{d\Phi(r)}{dr},$ (9)

where $\rho(r)$ is the cluster density profile. Using the equations of the
Universal Pressure Profile (see Section 2.1), we derived the following
expression for the dark matter potential:

$\Phi(r)=\frac{kT}{\mu
m_{u}}\ln\left[x^{\gamma}\left(1+\left(c_{500}x\right)^{\alpha}\right)^{\left(\beta-\gamma\right)}\right],$
(10)

this can be compared with dark matter potentials for the $\beta$-profile
(Krause, 2005) or the NFW-profile (Binney & Tremaine, 2008). The RMHD module
of PLUTO requires the input of the gravitational field strength to hold the
gas in place, this is derived using $\vec{g}=-\nabla\Phi$ to give:

$\vec{g}=-\frac{k_{B}T}{\mu
m_{u}c^{2}}\left[\frac{\left(\beta-\gamma\right)}{\left(R_{500}/(c_{500}r)\right)^{\alpha}+1}+\gamma\right]\frac{\vec{r}}{r^{2}}$
(11)

where $r=\sqrt{{x}^{2}+{y}^{2}+{z}^{2}}$ and $x,y$ and $z$ are distances in
the three Cartesian coordinate directions; $m_{u}$ is the atomic mass unit;
$c$ is the speed of light; $k_{B}$ is the Boltzmann constant and $T$ is the
temperature in Kelvin. In this model it is assumed that the gas has no self-
gravity and that it is held in place entirely by the dark matter potential.
The gas fractions ($f_{\mathrm{gas}}=M_{\mathrm{gas}}/M_{\mathrm{tot}}$) of
our model clusters range from $0.053-0.051$ (from low to high-mass) which are
at the lower end of the range of observed gas fractions (e.g. Wicker et al.,
2022) These low values may be in-part due to the debiased values of the UPP
parameters we employ which indicate gas pressure values (and so densities)
significantly lower than that predicted by the assumption of hydrostatic
equilibrium (see Section 2.2); they may also be due to a possible trend in gas
fraction with redshift, whereby lower redshift clusters have a lower gas
fraction and we have placed ours at $z=0$ (such a relation is hinted at by
Wicker et al. (2022)).

The nature of the UPP is that there is a cusp at the centre of the atmosphere;
Arnaud et al. (2010) highlight that at small radii, this is not realistic, and
that a reliable distribution is only expected from tens of kpc from the
centre. The difficulty in terms of coding this aspect of the atmosphere is
that the large gradient in pressure towards the centre, when discretised,
results in irregularities between the interpolated pressure distribution and
the body force implemented within PLUTO; these irregularities manifest
themselves as acceleration of the material, and so the atmosphere is not
balanced in the very centre. The way to get around this problem is to cut off
this unstable cusp at a suitable distance from the centre (we used a distance
of 4.2 kpc, determined by the resolution used) and created a zone of constant
value in that vicinity (i.e. set the pressure value within a certain radius to
be equal to the pressure value at that radius). The corresponding alteration
to the body force also needs to be made to stabilise this area (i.e. set the
gravitational field strength to zero). It should be noted that the injection
cylinder injects the jet at a distance of 6.3 kpc from the centre, and so the
region cut off does not impact the progression of the jet. We checked the
stability of our atmospheres by running our model, without a jet, for 200
simulation time steps (equivalent to $68.5\,$Myr); in this time a typical RMHD
jet would have progressed to the edge of the datacube. The flat distribution
at the very centre of the cluster creates a small discontinuity at its edge
which settles over time, otherwise the profile is static.

###  3.3 Simulation set up

For consistency, we use the code units of Paper 3. The simulation unit for
density, length and pressure are set to
$\rho_{0}=$3.01\text{\times}{10}^{-23}$\,$kgm-3; $l_{0}=2.1$ kpc and
$p_{0}=\rho_{0}c^{2}=$2.7\text{\times}{10}^{-6}$$ Pa, respectively. Finally,
the simulation unit for the magnetic field is calculated from
$B_{0}=c\sqrt{4\pi\rho_{0}}$, giving $1.84$$\,\mu$T.

The simulations are carried out on a static three-dimensional Cartesian grid
centred on the origin and extending to a length of $300$ kpc in each
direction. The central patch is a 4.2 kpc cube in width and is represented by
50 grid points in the y and z-directions and 10 grid points in the
x-direction; this provides sufficient resolution for the end of the injection
cylinder (see below). Either side of the central patch is a geometrically
stretched grid of 200 cells in the y and z-directions and 300 cells along the
x-direction. The resolution along the y and z-directions ranges from 0.084 kpc
at the centre to 6.9 kpc at the grid boundary; along the x-direction the
resolution ranges from 0.42 kpc at the centre to 2.1 kpc at the grid boundary.
The cell count is, therefore, $(n_{x},n_{y},n_{z})=(610,450,450)$; all outer
boundaries are set to ‘periodic’.

An injection cylinder is positioned in the centre of the grid, with the two
jets running along both directions of the x-axis. The advantage of this set-up
is that AGNs are bipolar outflows and it may be that the out-flowing jet from
one side will influence the jet on the other through back-flow of the radio
lobes (e.g. Krause & Camenzind, 2003; Krause, 2005; Antonuccio-Delogu & Silk,
2010; Hardcastle & Krause, 2013, 2014; Cielo et al., 2014). We use a jet
radius of 0.2 simulation units (0.42 kpc); a comparatively small value
compared to Paper 2, Paper 3 and Paper 4. As well as providing a more
realistic radius, this choice results in higher density injection for the same
jet power and this helps get the jet onto the grid and reduces the clouds of
ejecta surrounding the injection region (a problem identified in Paper 3);
this is particularly important for our current work as the UPP has a higher
density at the core than the equivalent $\beta$-atmosphere. We keep the
injection cylinder half-length the same as our previous work at 3.0 simulation
units (6.3 kpc).

Paper 4 highlighted the problem of material falling into the side of the
injection cylinder: this problem increases with the UPP as a result of higher
pressure values around the core. In order to solve this problem, we created a
thin layer around the curved surface of the cylinder (an annular cylinder) of
0.05 simulation units (0.105 kpc) and implemented a zero-gradient boundary
condition between it and the surrounding material for the magnetic field,
density, tracer and (critically) pressure (similar to the method employed by
Mukherjee et al. (2018, 2020)).

We used PLUTO version 4.4-patch2 for this study (Mignone et al., 2007); all of
the runs were carried out on the University of Hertfordshire High Performance
Computing facility. Each job was run on 384 Xeon-based cores, taking between
one and four weeks each. An output file was written by PLUTO every 50
simulation time steps (every 0.34 Myr).

We use the special relativistic magnetohydrodynamics (RMHD) physics module,
HLLD approximate Riemann solvers and a second order dimensionally unsplit
Runge-Kutta time-stepping algorithm, with a Courant-Freidrichs-Lewy number of
0.3. A divergence cleaning algorithm is used to enforce
$\nabla\cdot\textbf{B}=0$; we calculated the absolute value of the relative
divergence error (as defined by Pakmor & Springel (2013)) and found typical
values of $\sim${10}^{-3}$$. The model assumes a single-species relativistic
perfect fluid (the Synge gas) which is approximated by the Taub-Mathew
Equation of State (Taub, 1948; Mathews, 1971); for numerical stability
reasons, shock flattening was enabled through the use of a more diffusive
Riemann solver (HLL) and limiter (MIN-MOD); our simulations are non-radiative
for both jet and cluster material.

###  3.4 Model jet parameters

The jet is injected with a constant velocity of $0.994985c$; this corresponds
to a Lorentz factor of $\gamma=10$, well within the range observed in jets
(see Section 1). The jet power $Q_{\text{\scriptsize{RMHD}}}$ of a jet in SI
units, from English et al. (2016);

$Q_{\text{\scriptsize{RMHD}}}=\pi
r_{j}^{2}\nu_{j}\left[\gamma(\gamma-1)\rho_{j}c^{2}+\frac{\Gamma}{\Gamma-1}\gamma^{2}P_{j}+\gamma^{2}\frac{B^{2}}{2\mu_{0}}\right],$
(12)

where $r_{j}$ is the radius of the jet, $\rho_{j}$ is the jet density, $B$ the
magnetic field strength, $c$ is the speed of light, $\mu_{0}$ is the
permeability of free space, $\nu_{j}$ is the velocity of the jet (in units of
the speed of light), $\gamma$ is the Lorentz factor, $\Gamma$ is the adiabatic
index and $P_{j}$ is the pressure of the jet.

As in Paper 3 we limit our study to jets with equal contributions of enthalpy
and kinetic energy (i.e. the first and second terms of Equation 12). We inject
a helical magnetic field. The longitudinal component is set to a constant
value $B_{l}$ (for a particular jet power) and the toroidal component is
implemented as

$\text{B}_{y}=\text{B}_{t}(\text{z}/\text{r})\qquad\text{and}\qquad\text{B}_{z}=\text{B}_{t}(\text{y}/\text{r})$
(13)

for $r<r_{j}$ where $r_{j}$ is the radius of the jet. For this study, using
jet20_halo30 as the fiducial run, the jet values of $B_{l}$ and $B_{t}$ were
determined such that the values of the total magnetic energy of the lobes when
they have extended to an average length of $250$ kpc is $\sim${10}^{-2}$$ the
thermal energy (and $B_{t}$ and $B_{l}$ contribute equally). Values for other
runs were scaled up and down in proportion to the injected power, all values
can be seen in Table 3. Note that in this model the toroidal component
undergoes an unrealistic amplification at early times as the field lines are
stretched by the rapidly, outwardly expanding jet. This is suppressed by
gradually increasing the toroidal component from zero to maximum in the first
20 time steps (6.85Myr); from that point the injected field is constant. The
magnetic field evolution for the fiducial run can be seen in Figure 2. The
ratio of magnetic energy density to thermal energy density in the lobe for all
runs can be seen in Figure 7 where it can be seen that this value varies with
lobe length, jet power and mass of atmosphere; values range from
$\sim$2\text{\times}{10}^{-3}$$ to $\sim$3\text{\times}{10}^{-2}$$ for evolved
lobes across the range of parameters studied.

Table 3: Parameter values for the range of jet powers used in this study; power is measured in watts whereas the jet density ($\rho_{j}$), pressure ($P_{j}$), longitudinal ($B_{l}$) and toroidal ($B_{t}$) magnetic fields are all measured in simulation units. Each power is run into an atmosphere (marked as YY here) of either $03,10,30$ or $90$ which represent atmospheres of $M_{500}=0.33,1,3$ and $9\times${10}^{14}$\,\text{h}^{-1}_{70}M_{\odot}$ (see Table 2). run | power | $\rho_{j}$ | $P_{j}$ | $B_{l}$ | $B_{t}$
---|---|---|---|---|---
|  ${10}^{38}$W |  ${10}^{-6}$sim. |  ${10}^{-6}$sim. |  ${10}^{-3}$sim. |  ${10}^{-3}$sim.
jet05_haloYY | $0.5$ | 0.6452 | 0.2043 | 0.280 | 1.350
jet10_haloYY | $1$ | 1.2903 | 0.4091 | 0.403 | 1.909
jet20_haloYY | $2$ | 2.5806 | 0.8175 | 0.570 | 2.700
jet40_haloYY | $4$ | 5.1612 | 1.6358 | 0.806 | 3.818
Figure 2: Magnetic field evolution in the lobe (LHS) and the ratio between the
magnetic and thermal energy in the lobe (RHS) for the fiducial run
(jet20_halo30). The magnetic field of the jet is injected with both toroidal
and longitudinal components (i.e. helical injected field); the values were
chosen so that the total magnetic energy of the lobes is $\sim${10}^{-2}$$ the
thermal energy, and the toroidal and longitudinal components are roughly equal
in terms of magnetic field energy, once the lobes have reached an average
length of $250$ kpc (see text).

###  3.5 Postprocessing

Following Paper 3, the jet is injected with a conserved tracer quantity of
value of 1.0 (and zero elsewhere). Lobes are defined by tracer values
$>10^{-3}$. The bow shock surface (between the shock and the undisturbed
ambient medium) is identified in a similar way to the tracer, by line-tracing
from the edges of both sides of the volume towards the centre and finding
where the radial velocity exceeds the defined value of $75$ kms-1.

##  4 RESULTS

###  4.1 Dynamics

Figure 3 shows density images for all runs when they have progressed to 250
kpc, these are taken as slices through the centre of the xy-plane.
Qualitatively, we observe several features also present in our previous
simulations, particularly those of Paper 3. We observe that some simulations
have sufficient buoyancy to lift lobe material clear of the central regions of
the simulation, in particular the higher-mass atmospheres and the lower-power
jets. The higher mass atmospheres have greater values of gravitational
potential and so will exert a greater force on the rising bubbles - but the
lower power jet simulations can clear the central regions of lobe material
simply because they take longer to reach 250 kpc and so the buoyancy forces
act for longer. It is clear that the axial ratio varies significantly across
these simulations; higher mass atmospheres and lower power jets create lower
values for axial ratio by inflating less elongated prolate spheroidal bubbles
than higher power and lower mass atmospheres. Many of our runs have asymmetric
lobes, this is due to the slight density perturbations in the initial
environment of our models; these are more pronounced than in our previous work
as a result of the lower power jets used here (lighter jets are more
susceptible to their path being altered by such inhomogeneities). Lobe
asymmetries are a feature of real lobes (e.g. Hardcastle et al., 1997b).

Figure 3: Logarithmic density images (xy-slice) for the full range of runs,
each presented once the simulated radio galaxy has extended to an average lobe
length of 250 kpc. Cluster mass $M_{500}$, jet power $P$ and time to reach an
average lobe length of $250$ kpc, $t_{250}$ are shown on each image.

Figure 4 shows the evolution of the lobe dimensions for all runs in this study
and has much in common with both the non-relativistic runs of Paper 2 and the
mildly-relativistic runs of Paper 3. The progression of the jet through the
cluster is governed by the balance between the lobe pressure and momentum flux
of the jet with the density and pressure of the ambient material at the end of
the lobe. As expected, higher power jets progress faster through the same
environment as a result of the greater density and so momentum flux of such
jets (given we have a fixed injection speed). Similarly, the richer
environments result in the same power jet progressing slower as a consequence
of the higher ambient pressure and density of the cluster material. The
evolution of lobe length is more linear than our previous work with $\beta$
models; this is because the $\beta$-profile has a near-uniform density at the
core which slows the jet until it reaches a steeper density slope further out;
whereas the UPP has a decreasing density slope throughout the run as a result
of the cusp at the centre. As expected, the higher power jets accelerate away
from the more dense core, but we also found that the higher mass atmospheres
(provided the power of the jet was low enough) resulted in the shock front
decelerating at longer lengths - corresponding to the inflated lobes seen in
the density images of Figure 3.

The bottom panel of Figure 4 demonstrates that these simulations are not self-
similar - a constant value for volume/length3 would be expected for self-
similarity. This finding agrees with our previous work, is a feature of other
simulations (e.g. Krause, 2005) and observations (e.g. Hardcastle & Worrall,
2000). It is also clearly seen in the density images of Figure 3 where low
power jets moving into rich environments progress more slowly than high power
jets moving into poor, and so have lobes of smaller aspect ratio (i.e. are
fatter). We would expect the slower jets to have a smaller aspect ratio as
they take longer to reach 250 kpc and so the transverse expansion (which is
near-adiabatic) progresses further. It seems likely that KH instabilities also
play a role in determining the aspect ratio as the lobes have time to lift
away from the cluster core and expose the fast-flowing jet to the near-
stationary ambient medium, the jet is no longer protected by the lobe and KH
instabilities set in. These instabilities disrupt the base of the jet and
effectively increase the opening angle of the jet-launch, resulting in a
larger working surface at the jet termination and so decreases the advance
speed of the lobe further. The impact of the disruption to the jet caused by
KH instabilities can be seen in Figure 17.

Figure 4: Evolution of lobe length (top); volume (middle) and the ratio of
volume to the cube of lobe length (bottom) for all runs.

###  4.2 Energetics

The total thermal, kinetic, magnetic and potential energies contained in the
lobe and shocked regions can be plotted against time to investigate the
efficiency of the energy injection process. An example is presented in Figure
5 which is for a jet of power = $1\text{\times}{10}^{38}$W running into an
atmosphere of $M_{500}=3\times${10}^{14}$\,\text{h}^{-1}_{70}M_{\odot}$. The
total energy measured in the simulation is almost identical to that injected
into the system, as expected.

Figure 5: Energies in the lobe and shocked regions as a function of time for
jet10_halo30.

The ratio of the energy stored in the shocked region to that stored in the
lobes (Figure 6) shows us that more powerful jets and lower mass atmospheres
both result in a greater fraction of the energy going into the shocked region.
The values seen here are generally higher than those of our previous studies;
this is likely to be a result of using what are effectively lower mass (and
more realistic) atmospheres than those used previously. The values in this
study (ranging from about $1$ to $4$) are double the values we observed in our
previous work, although smaller than those found in the simulations by Perucho
et al. (2011, 2014, 2019, 2022) who find values in excess of $10$; however,
their set-up is considerably different to ours and involves a jet perturbation
which may influence the transfer of energy from the lobe to the shocked
region. In Paper 1 and in particular Paper 2 we identified that the
environment played a role in determining the ratio of shocked to lobe
energies: in particular, the greater the density slope the greater the
fraction of energy found in the shocked region (i.e. for the $\beta$-profile
this was measured by the parameter $\beta$). For our UPP runs, the density
slope is much greater at the core of the cluster and so our current, higher,
values can be explained as following this same trend. As highlighted in
Section 4.1, KH instabilities impact the lower-power jets moving into richer
environments more and result in fatter lobes. Some studies suggest that such
turbulence in the ICM does not contribute significantly to the heating of the
ICM (e.g. Weinberger et al., 2017; Bourne & Sijacki, 2017; Martizzi et al.,
2019) whereas others say that they do (e.g. Gilkis & Soker, 2012; Hillel &
Soker, 2016, 2017). We conclude that the turbulence in our models does not
apper to significantly contribute towards heating of the ICM as our models
with the greatest jet disruption result in a smaller fraction of the total
energy going into the shocked region. An alternative explanation is provided
by Tang & Churazov (2017): in their simulations of rapid outbursts they found
that more explosive emission events resulted in a greater proportion of the
injected energy going into shocks, with up to $\sim 88\%$ for
instantaneous/explosive emission and $\sim 0\%$ for the longest and most
gentle injection. The rapidity of the injection, therefore, plays a major role
in determining the significance of the heating effect; this agrees with our
results as our high-power jets moving into poor environments have the most
rapid expansion and the greatest heating effect.

Figure 6: The ratio of the energy contained in the shock to that contained in
the lobe, against lobe length for all runs. Low power jets in the highest mass
atmospheres (jet05_halo90 and jet10_halo90) can be seen to have erratic lines
here; this is due to their low shock advance speed such that the method of
detection (speed falls below $$2.5\text{\times}{10}^{-4}$c$) is above the
advance speed in the outer cluster and so the shock front is not identified
accurately; full runs included here for completeness.

The ratio of magnetic energy density to thermal energy density (Figure 7)
varies with power, atmospheric mass and lobe length in a systematic way: all
lines appear to be roughly cubic in form; lower power jets have greater values
at shorter lobe lengths and then decrease with increasing lobe length, whereas
higher power jets increase monotonically. Atmospheric mass alters the ratio
such that for all powers, higher mass increases the ratio at shorter lengths,
but decreases it at longer; there is a crossover point on the graph where all
atmospheres have the same ratio (roughly) for a particular power. Through
their work in comparing models with observations, Böhringer et al. (2016) find
that the magnetic energy density is two orders of magnitude smaller than the
thermal energy with the ratio falling in the range
$5-$7.5\text{\times}{10}^{-3}$$; this means that the magnetic fields in
clusters are considered not to be dynamically important (i.e. plasma beta
$\equiv P_{\text{th}}/P_{\text{B}}\sim 100$). Figure 7 indicates that the
magnetic field is not dynamically important, given that all runs have values
$\sim 0.01$, although it is possible that it could be locally dynamically
important.

Figure 7: The ratio of magnetic energy density to thermal energy density
variation with lobe length for all runs.

###  4.3 Synchrotron emission and polarimetry

We simulate Stokes $I$, $Q$ and $U$; and from these we calculate the
polarisation and the magnetic field direction. We follow the method employed
in Paper 2 and Paper 3; a summary of this can be found in Appendix B.

By integrating the synchrotron emission over the whole of the source for each
simulation data cube we produce graphs for the total luminosity against lobe
length; Figure 8 shows the results of such an analysis for the full range of
power and atmospheric mass used in this study when viewed at various angles to
the jet axis. In line with similar plots from Paper 2 and Paper 3, at
$90^{\circ}$ the luminosity initially increases with lobe length and then
begins to level out, sometimes reaching a maximum or a plateau. In keeping
with Paper 2 we find that that atmosphere used impacts the luminosity (in
Paper 2 the parameters $r_{c}$ and $\beta$ for the $\beta$-profile were
varied). The general trend is that for higher mass atmospheres the higher the
synchrotron output. On our $90^{\circ}$ chart we have included horizontal
lines based upon the relation between jet power and the synchrotron luminosity
at 151 MHz by Willott et al. (1999); this relation predicts that $Q\propto
L_{151}^{6/7}$. Here we have normalised this to our results for jet power
$1\text{\times}{10}^{38}$ W (ignoring the highest mass) at late times. Whilst
there is reasonable agreement, our results suggest a relation with a power
greater than 6/7. The fact that our results differ to those of Willott et al.
(1999) is to be expected as they assumed a power-law atmosphere in their model
whereas we used the UPP; and their relationship does not capture the time
evolution that we observe. Clearly this relation would also have to be
adjusted to incorporate viewing angle if it were to be used on all our plots.
Varying the angle of view has a significant impact on the emission luminosity,
particularly for angles $\sim 1/\gamma$ and less (where $\gamma$ is the
Lorentz factor) as Doppler boosting of the jet becomes significant (about
$6^{\circ}$ in our study) and the emission from the jet dominates over that of
the lobe. We must remind ourselves that we have made a number of assumptions
about the pressure, magnetic field structure and electron energy spectrum of a
real jet, and so the results we present here for viewing the jet at small
angles must be regarded with caution; nevertheless, end-on jets are observed
(e.g. blazars, Blandford et al. (2019)) although we feel that the amount of
Doppler boosting is exaggerated in our models. Movies showing the simulated
synchrotron emission at various angles and lobe lengths can be found at
http://uhhpc.herts.ac.uk/~ms/synchrotron.html.

Figure 8: Evolution of synchrotron luminosity with lobe length for all runs
viewed at various angles to the jet axis as indicated: 90∘ (top left), 30∘,
20∘, 10∘, 5∘ and along the jet axis (0∘). The black horizontal dashed lines on
the top-left plot represent the emission predicted by Willott et al. (1999),
which has been normalised to the $1\text{\times}{10}^{38}$W lower cluster mass
lines at late times (see text). The line for lowest power run in the highest
mass atmosphere (jet05_halo90) becomes erratic once it has progressed beyond
150 kpc; this jet finds difficulty in punching its way through the ambient
medium and the tracer mixes more than other runs and falls below the threshold
of ${10}^{-3}$, and so the location of the lobe becomes unreliable and is
included here for completeness. Figure 9: Sequence of synchrotron emission
maps; top image observing down the jet axis, each subsequent image is
positioned an additional $10^{\circ}$ from the axis, the bottom image is
looking down the counterjet. Images for a jet of power
$2\text{\times}{10}^{38}$W running into an atmosphere of
$M_{500}=9\times${10}^{14}$\,\text{h}^{-1}_{70}M_{\odot}$ with jet extended to
an average length of 250 kpc; arbitrary logarithmic brightness scale. Doppler
boosting of the jet and counterjet is visible at small angles to the jet-
counterjet axis; note that the logarithmic scale makes this appear less
pronounced; the true magnitude of the Doppler boosting of the jet is more
easily seen in Figure 8.

We created simulated synchrotron two-dimensional emission maps; an example can
be seen in Figure 9 where the angle of view has been progressively increased
in order to observe the simulated Doppler boosting. We also created emission
maps for Stokes $Q$ and $U$ and $P$ (polarised intensity); examples can be
seen in the top panel of Figure 10 and are very similar to those produced by
our previous work in Paper 2 and Paper 3, although a greater level of fine
detail is revealed as a consequence of the higher resolution used in the
current study. We can overlay a synchrotron emission map with vectors with
magnitude that represents the magnitude of the fractional polarisation and
with a direction corresponding to that of the magnetic field. Such an image
viewed at $90^{\circ}$ to the jet axis is presented in the middle panel of
Figure 10 and shows the same trends found in to our previous work in that the
evolved lobe polarisation is higher at the edges where the magnetic field
direction is parallel to the lobe-shock boundary. Furthermore, in the bottom
panel of Figure 10 we demonstrate that these findings are also true when
observing from an angle (here at $30^{\circ}$ to the jet axis), in common with
what was found in the simulations of Huarte-Espinosa et al. (2011) and seen in
observations by Hardcastle et al. (1997b). These similarities give us
confidence that the conclusions drawn from our previous work remain valid for
these more relativistic and higher resolution studies which use a more
realistic helical magnetic jet injection.

Figure 10: Synchrotron emission maps for jet20_halo90 with lobe extended to
250 kpc. Top panel: (from top left) Stokes $I$ (synchrotron emission), $Q$,
$U$ and the polarised intensity $P=\sqrt{Q^{2}+U^{2}}$ (bottom right). Middle
and Lower panels show synchrotron emission maps overlaid with vectors
representing the direction of the magnetic field and the magnitude of the
linear polarisation. Middle panel viewed at $90^{\circ}$ to the jet axis and
Lower panel at $30^{\circ}$. Arbitrary brightness scale used.

##  5 DISCUSSION

###  5.1 Comparisons with the $\beta$-model

Here we draw comparisons between our UPP models and $\beta$-models in order to
establish whether our UPP model leads to different results from the
established $\beta$-profile. We also compare our results with observations of
Doppler boosting and observations of radio jets, and report on our resolution
study.

There is no established way to create an equivalent $\beta$-profile from a UPP
profile: Hardcastle (2018) present a graph comparing these two profiles, and
it is clear that they differ considerably at the cluster centre in particular.
Starting with a UPP profile, we derived the equivalent $\beta$-profile using
an iterative technique beginning with the general profile of $\rho_{0}=1.0$
simulation units, $r_{c}=0.1R_{500}$ and $\beta=0.5$; and then varied each of
these three parameters by a small amount and minimised the difference between
the two profiles. The UPP profile has a ‘core radius’, defined by Arnaud et
al. (2010) as $r=0.1R_{500}$; there is no indication that this is equivalent
to the $\beta$-profile $r_{c}$ but it gives an order of magnitude value to
start with; the initial $\rho_{0}$ and $\beta$ are typical values employed by
Paper 2 and Paper 3. The ‘best fit’ $\beta$-model values we derived for each
UPP model atmosphere can be seen in Table 2.

Density profiles of the UPP and corresponding $\beta$-model for the range of
cluster mass used in this study can be seen in Figure 11. The two models
differ most at the cluster centre and follow very similar curves out to the
limit of this study (i.e. 300 kpc). Beyond this range of study, the cluster
density is comparatively low for both models and so any variation at large
distances would have an insignificant impact on the progression of the jet and
lobe in comparison with moving through the higher density regions closer to
the core.

Figure 11: Fitting the $\beta$-profile to the UPP by minimising the difference
between the two curves (by varying the $\beta$-model parameters $\rho_{0}$,
$r_{c}$ and $\beta$). Continuous lines are for the UPP atmospheres and dashed
lines are for the corresponding $\beta$-profile. Note that the log-log plot
used here highlights the difference at the core.

Govoni et al. (2010) provide $\beta$-model parameters for the much-studied
Coma cluster of $\rho_{0}=$3.5\text{\times}{10}^{-3}$\,$m-3, $r_{c}=245\,$kpc
and $\beta=0.654$. We would not expect these to be exactly replicated in an
‘average model’ such as the UPP, although one would hope that values would be
the same order of magnitude. Assuming a Coma cluster mass with
$M_{500}=$2.9\text{\times}{10}^{14}$\text{h}^{-1}_{70}M_{\odot}$ (Lyskova et
al., 2019), and so using our cluster mass of
$M_{500}=$3\text{\times}{10}^{14}$\text{h}^{-1}_{70}M_{\odot}$, as well as the
general temperature profile used in this study (see Section 2.3) which has a
central temperature value of 2.4 keV, we read off the parameters from Table 2
and converting from simulation units we obtain:
$\rho_{0}=$1.4\text{\times}{10}^{-3}$\,$m-3, $r_{c}=95\,$kpc and $\beta=0.54$
and so we do have the correct order of magnitude and values derived from the
UPP which differ from observed values by no more than a factor of $\sim\\!3$.

We ran a jet of power $1\text{\times}{10}^{38}$W into a
$M_{500}=3\times${10}^{14}$\,\text{h}^{-1}_{70}M_{\odot}$ atmosphere
(jet10_halo30) and the same power jet into the equivalent $\beta$-atmosphere
(beta10_30) (see Table 2 for parameter values). Both runs extended to 250 kpc
and density images can be seen in Figure 12 with dynamic and energetic
comparisons presented in Figures 13 and 14. The atmosphere model chosen
impacts the morphology of the lobe as well as its dynamics and energetics. The
UPP is more able to clear lobe material from the central regions and so leads
to a more distinct pair of bubbles whereas the $\beta$-model forms one
elongated shape where the lobes remain merged at the centre. However, given
the rich variety of observed lobe shapes, examples can be found which
correspond to both shapes seen here: 3C 438 has lobes which are more separated
whereas 3C 20 has lobes which merge near the core (radio images of these lobes
and others can be seen in Hardcastle et al. (1997b)). The $\beta$-model has a
higher aspect ratio at early times (Figure 13) although differences become
negligible once lobe lengths of hundreds of kpc are reached; this is
consistent with the atmosphere models given that the major difference between
$\beta$ and UPP is at the core. Comparisons are more marked for the evolved
lobes when considering the energetics: the top chart of Figure 14 shows that
the lobe in the $\beta$-model passes on energy to the shocked region less
efficiently at early times but ends up transferring the greatest by 250 kpc.
Similarly the magnetic/thermal ratio of the lobe is lower at early times but
ends up being greater at late times (Figure 14 (bottom)). The values for
individual lobes have been plotted on Figure 14 (and 13) and these give an
indication of the variability within these runs; this demonstrates that there
is a genuine difference between the UPP and $\beta$-model atmospheres which
impacts their energetics, particularly at late times.

Figure 12: Logarithmic density images (xy-slice) for a jet of power
$1\text{\times}{10}^{38}$W and moving into a UPP atmosphere of
$M_{500}=3\times${10}^{14}$\,\text{h}^{-1}_{70}M_{\odot}$ (top) and the same
power jet moving in to the equivalent $\beta$-atmosphere (bottom). Parameters
used to create these two atmospheres can be seen in Table 2. Both models shown
have progressed to 250 kpc and take almost exactly the same time to do this
(56 Myr for jet10_halo30 and 57 Myr for beta10_halo30). Figure 13: Evolution
of length (top) and the ratio of volume to length3 (bottom) for jet10_halo30
(UPP run), beta10_halo30 ($\beta$ run); jet10_halo30_highres (UPP high
resolution run) and jet10_halo30_lowres (UPP low resolution run). See Figure
12 for density images of jet10_halo30 and beta10_halo30. The thick line is the
average for the two lobes and the thin lines are for individual lobes; this
gives a sense of the scatter between the lobes. Figure 14: Evolution of the
ratio of shocked region to lobe region energy (top) and lobe magnetic to
thermal energy (bottom) for jet10_halo30 (UPP run), beta10_halo30 ($\beta$
run); jet10_halo30_highres (UPP high resolution run) and jet10_halo30_lowres
(UPP low resolution run). See Figure 12 for density images of jet10_halo30 and
beta10_halo30. The thick line is the average for the two lobes and the thin
lines are for individual lobes; this gives a sense of the scatter between the
lobes.

###  5.2 Synchrotron imaging

As seen in Section 4.3, our model reproduces the expected Doppler boosting of
the jet when viewed along (or at small angles to) the jet axis. This effect
can be measured by the brightness ratio ($I_{j}/I_{cj}$) and a plot of jet
base sidedness can be created from observations by rotating the synchrotron
image of a bipolar jet by $180^{\circ}$ about the nucleus and dividing the jet
image by the counterjet (e.g. Laing, 1996; Hardcastle et al., 1997a). We copy
this methodology for the simulated emission of jet05_30 viewed at $30^{\circ}$
to the jet axis to obtain Figure 15 (top). This ratio of jet base sidedness
can also be derived from the jet bulk velocity $\beta$ and viewing angle
$\theta$ as (e.g. Laing, 1996)

$\frac{I_{j}}{I_{cj}}=\left[\frac{1+\beta\cos\theta}{1-\beta\cos\theta}\right]^{2+\alpha}$
(14)

where $\alpha$ is the synchrotron spectral index (taken to be 0.5). This
formula assumes isotropic emission in the rest frame and so can only be used
as an approximation. From our simulation we calculate the brightness ratio
using both the left hand (simulated synchrotron brightness ratio) and right
hand (bulk velocity) sides of Equation 14 to produce Figure 15 (bottom). The
line derived from the bulk velocity falls from the initial maximum more
smoothly than that derived from the brightness ratio, this is to be expected
as an isotropic magnetic field has been assumed for the bulk velocity
calculation, whereas the brightness ratio calculation incorporates the
variation in the simulated magnetic field. A very similar evolution for the
simulated brightness ratio was found by Laing (1996) (their Figure 4) from
their synchrotron observations of 3C 31, a decelerating jet. Measurements of
the bulk velocity of our jet show that it is also decelerating in this inner
region; this is as a result of buoyancy having lifted the lobe and exposed the
jet to the ambient material which is entrained into the edge of the jet and so
slowing it down.

Figure 15: A grey-scale image of jet sidedness (top) for jet05_halo30 viewed
at $30^{\circ}$ to the jet axis. It is constructed by rotating the simulated
synchrotron image by $180^{\circ}$ about the nucleus and dividing by itself.
Values above 5 are saturated in this image. Longitudinal jet sidedness
profiles (bottom) derived from the simulated image of jet sidedness (shown
above) and also from the simulated bulk velocity measures (see text).

Figure 16 shows an image of the radio galaxy 3C 349 (reproduced from
Hardcastle et al. (1997b)); this can be compared with the synchrotron image of
Figure 10 (top left), the images of Figure 9 or the density images of Figure
3. Favourable similarities can be noted including: the approximate axial
ratio, the effect of buoyancy which thins out the radio-emitting lobe from the
central regions, the rough large-scale symmetry and the overall length. Note
that the example chosen (3C 349) is imaged close to $90^{\circ}$ to the jet
axis as otherwise Doppler boosting would complicate comparisons with the
density images presented here. The presence of Kelvin-Helmholtz instabilities
on the lobe surface and the small-scale morphological variability is evident
in our simulated density and synchrotron images as a result of density
inhomogeneities in the ambient cluster; these are also evident in the observed
radio image. In particular, some of our lower power runs produce fatter lobes
in the vicinity of the cluster centre, but then the lobe growth increases once
it reaches lower density material in the outer cluster and subsequently
produces a ‘neck’, a similar feature is seen on the upper lobe of 3C 349.

Figure 16: 3C349 at 1.80-arcsec resolution. Contours at
0.25$\times$(-2,-1,1,2,4,….) mJy beam-1. Image taken from Hardcastle et al.
(1997b).

###  5.3 Lack of hotspots

As discussed in our earlier work (Paper 3), our models do not reproduce the
hot spots characteristic of ‘edge brightened’ FR II radio galaxies (these are
clearly visible in observations such as Figure 16). Knots and hotspots are
synchrotron emission created by highly energetic particles; these features are
normally far from the injection and recollimation sites near the central
nucleus such that adiabatic losses (Longair et al., 1973) and spectral ageing
(e.g. Kardashev, 1962) would prevent particles of sufficient energy reaching
them so that there needs to be local acceleration of particles, such as the
process of diffusive shock acceleration (DSA, Bell (1978)). We mapped the
internal Mach number of our simulations: $\Gamma\beta/\Gamma_{s}\beta_{s}$;
where $\Gamma$ is the Lorentz factor, $\beta$ is the bulk velocity in natural
units and the subscript $s$ denotes values for the sound speed (see Figure
17). We found that the jet has a Mach number in excess of 1 for almost all its
length, and so a shock could take place along the jet right up to the tip of
the lobe. Simulations by Massaglia et al. (1992, 1996, 1997) indicate that the
evolution of Kelvin-Helmholtz (KH) instabilities at the surface of the jet
results in the formation of shocks which may be the cause of knots. Such
disruption of the jet through KH instabilities may also be a factor
influencing the termination of the jet in wide-angle tail radio sources
(Hardcastle & Sakelliou, 2004; Hardcastle et al., 2005), particularly where a
hotspot is seen at the base of the plume. In our simulations we can see that
shortly before the end of the lobe the Mach number falls to subsonic levels as
the jet encounters decelerating material back-flowing into the lobes; and the
jet terminus moves around in these snapshot images, in keeping with the
description of Scheuer (1974) and his ‘dentist drill’ model; we see that this
feature is increased when a jet of the same power is injected into a richer
environment as this reduces the Mach number and so leads to a less stable jet.
We note also that in the richer environments the lobe lifts away from the core
of the AGN through buoyancy and exposes the outer surface of the jet to the
higher-density ambient material; this should reduce the stability of the jet
and in our models we see significant disruption through entrainment and
deceleration of the outer surface of the jet leading to a ‘sheath’ of slower-
moving material around the ‘spine’ of the jet; the slower moving sheath is
more susceptible to KH instabilities. We suspect that the creation of this
sheath of material may be due to our boundary conditions when launching the
jet, resulting in an amplification of the KH instabilities beyond what would
be realistic.

Figure 17: Internal Mach number ($\Gamma\beta/\Gamma_{s}\beta_{s}$) for a jet
of power $0.5\text{\times}{10}^{38}$ W running into an atmosphere of $M_{500}$
(measured in units of $\times${10}^{14}$\,\text{h}^{-1}_{70}M_{\odot}$): $0.3$
(top left), $1$ (top right), $3$ (bottom left) and $9$ (bottom right). Images
are of x-y slices through the midpoint and the plots are of the maximum Mach
number in the y-z plane at that distance from the core (to take into account
the movement of the jet out of the x-y slice shown).

###  5.4 Aspect ratio

The aspect ratio of our models varies significantly (see Figure 3). FR II’s
tend to form prolate spheroidal-shaped lobes with an axial ratio between $1.3$
and $6$ (Leahy & Williams, 1984) with the largest fraction between 1.5 and 2
(Mullin et al., 2008), where axial ratio is defined as the ratio of length to
width. Our lower power jets moving into richer environments fall well within
these ranges, although our higher-power jets are at and beyond the upper
range. The jet power and cluster masses of our models match well those seen in
observations, and so there must be another feature of high power jets which is
missing from our models. Numerical simulations limit the resolution by means
of a minimum grid size and so reduce the growth of instabilities experienced
by the model jet as it moves through the material within the lobe which would
otherwise slow the jet’s progress (e.g. Mignone et al., 2010); this effect
explains why the lobe growth is often slowed when higher resolution
simulations are conducted in convergence testing (e.g. Krause & Camenzind,
2001; Perucho et al., 2019). Sometimes a method is employed to compensate for
the reduced instability growth by limiting the flow of the jet so the model
has time to form realistic-shaped lower aspect-ratio lobes. One such technique
is to introduce a helical perturbation to the injection-velocity; this is
achieved through applying small velocities normal to the jet flow at the
injection cylinder, for example in the simulations of Perucho et al. (2019,
2022) where a sum of normal sinusoidal velocity perturbations of different
angular frequencies is used. They observe a very stable jet (phase I) up to
100 kpc, and beyond this (phase II) a ‘dentist drill’ phase whereby the jet is
appreciably disrupted; they refer to this second phase happening in their
model as a result of the ‘wobble’ introduced and cite irregularities near the
jet head as evidence of helical oscillations (e.g. Harwood et al., 2016,
2017). The phase II results in a decrease in the advance speed and subsequent
decrease in aspect ratio. In our simulations we do not employ such a model,
instead relying upon the numerical noise of our simulation and the
asymmetrical fluctuations in density created at set-up to deviate the jet as
it traverses the non-uniform atmosphere; whilst this is sufficient for low-
power jets, it was not sufficient to prevent the higher power runs from
forming unrealistically high aspect ratio lobes.

However, it is well known that at least some powerful FR II radio sources are
at the same time large, have a stable jet, and a fat lobe (e.g., Pictor A,
Hardcastle et al. (2016)). A straightforward explanation of the fat lobes may
then be jet precession. Other indicators of jet precession have been
identified in recent simulations (Horton et al., 2020a; Horton et al., 2020b;
Giri et al., 2022) and are frequently seen in powerful radio sources (Krause
et al., 2019). It appears difficult to explain these findings otherwise. While
it has been widely demonstrated that lighter jets inflate fatter lobes (e.g.
Krause, 2005; Gaibler et al., 2009; Guo, 2015, 2016; Su et al., 2021;
Weinberger et al., 2023) and that interaction with the an inhomogeneous ICM
can significantly alter the morphology of the lobe (e.g. Wagner & Bicknell,
2011; Wagner et al., 2012; Mukherjee et al., 2018, 2021a), lowering the jet
density at constant power would require an increase in the Lorentz factor. We
already use a Lorentz factor of ten, which is towards the upper end of
plausible values, as constrained, e.g., by multi-epoch Very Long Baseline
Interferometry data. Hence our jet densities could not plausibly be much
lower. We are therefore left with jet precession as the best explanation for
fat radio lobes in galaxy clusters.

###  5.5 Resolution study

An increase in resolution tends to result in a simulation which takes us
closer to physical behaviour and closer to convergence (e.g. Krause &
Camenzind, 2003; Perucho et al., 2019). In their resolution study of 3D RMHD
jets Mignone et al. (2010) found that higher resolution jets were able to
maintain their Lorentz factor better than equivalent lower-resolution runs;
whereas Weinberger et al. (2023) showed that the identification of convergence
of their models was non-trivial as the factors influencing the length of the
jet depended upon the resolution used, with lower resolutions being influenced
more by numerical diffusion and higher resolution more by KH instabilities
resulting in a turnaround in their jet distance versus time plot, also compare
with Krause & Camenzind (2001). Through the use of stretched grids we have a
higher resolution at the centre of our simulation than our previous studies:
In comparison with the uniform grids of ;  and Paper 4 the resolution in the
x-direction is 3.75 times higher and in the yz-plane is 18.75 times higher at
the centre. In order to judge to what degree resolution influences our models
we ran simulations with different numbers of grid points in the x-direction.
In our normal-resolution runs, we have a geometrically stretched zone of 300
cells on either side of the central uniform zone; in our low resolution run we
replaced this with a zone of only 150 cells on either side of the central
uniform zone, which essentially halved the number of grid points in the
x-direction. In our high resolution run we increased the stretched grid to 300
cells on either side of the central uniform zone, and so essentially doubled
the resolution along the x-direction. We compared the progression of the lobe
and the lobe volume/length3 ratio for jet10_halo30 for all three resolutions
(see Figure 13) and found very similar results, particularly for the two
highest resolutions. Whilst not identical, the fact that the runs all cross in
a number of places suggests that the variability between the runs is slight
and may be influenced not so much by resolution but by the interpolation of
mapping the initial conditions onto a different grid at startup, impacting the
exact evolution of the turbulence. However, our energy analysis (Figure 14)
indicates a more significant difference: the lower resolution run follows a
separate path, whereas the higher resolution average lines often fall within
the variation indicated by both of their separate lobes, particularly away
from the cluster centre. We conclude from this limited-resolution study that
models of evolved lobes are very close to being converged in terms of both
dynamics and energetics with models of twice the resolution along the x-axis.

##  6 SUMMARY AND CONCLUSIONS

We created numerical models of realistic cluster atmospheres by utilising the
Universal Pressure Profile (UPP) of Arnaud et al. (2010), applied a typical
temperature profile of a cool core cluster, and ran relativistic helically-
magnetised jets into them. We modelled cluster atmospheres over the range of
typically-observed masses, from the group-cluster boundary up to high-mass
examples and our jet powers varied from the FR I/FR II boundary up to high-
energy values. In keeping with our previous work, this cluster model also
incorporates a multi-scale, tangled magnetic field which is related to the
cluster density profile (Huarte-Espinosa et al. (2011), Paper 2, Paper 3 and
Paper 4).

This study has built upon our previous work and improved some aspects of those
models, moving closer to realism. Our use of stretched grids has enabled a
much greater resolution along the jet axis and so a smaller and therefore more
realistic jet injection cross-sectional area (100x smaller than that used
previously). This has removed the unrealistic plumes of lobe material which
were created around the cluster core in our previous models (see Paper 3). In
the current work we have employed relativistic jets with a Lorentz factor of
$\gamma=10$, well within the established range of values observed for real
jets, whereas our previous relativistic studies (Paper 3 and Paper 4) used a
value of $\gamma=2.7$, which is very much at the lowest boundary of the range
of realistic values. In addition, our use of a helical magnetic field (rather
than injecting a purely toroidal field) is another step closer to realism.

In keeping with our previous work, we find that higher power jets inflate
lobes with a larger aspect-ratio, progress faster and create larger volume
lobes for the same age. In addition, we again demonstrate that the cluster
environment has a significant influence on the evolution of the lobe and we
see that atmospheres with a larger mass (measured by $M_{500}$) inhibit the
growth of the lobes more and result in a lower aspect ratio. In our previous
work we found that the ratio of shocked energy to lobe energy was a constant,
whereas here we find that this increases with jet power and decreases with
cluster mass. Our results for the total synchrotron luminosity follow a very
similar pattern to that of our previous work (increases with jet power)
although previously we found little variation with atmosphere; here we find
that lower power jets are impacted more by the cluster mass, the general trend
being that lower mass atmospheres result in greater luminosity. It was seen
that, within the range of parameters studied here, the jet power has a much
greater influence on the total synchrotron luminosity than cluster mass.

Synchrotron imaging, maps of the Stokes parameters and charts indicating the
direction of the magnetic field are very similar to those we created in our
previous work, showing the filamentary structure in a greater level of detail
as a result of the higher resolution in the centre of the model; favourable
comparisons can also be made with the simulations of other researchers (e.g.
Matthews & Scheuer, 1990; Tregillis et al., 2001) and with observations (e.g.
Saikia & Salter, 1988; Hardcastle et al., 1997b). In this current study we
have also seen the impact of Doppler boosting on the synchrotron image and
created sidedness plots which resemble those derived from observations of 3C
31.

Our simulations are the first to use the UPP in a numerical simulation of jet
feedback and the resultant density and synchrotron images of the lobes compare
very well with images of radio galaxies from observations. Using the same jet
parameters, we created an equivalent $\beta$-profile atmosphere and conclude
that whilst the dynamics of the UPP and $\beta$-profile are very similar
(particularly at large distances from the core), the greatest differences lie
in the energetics where significant variations can be seen in the distribution
of energy between the lobe and shock regions and also between the magnetic and
thermal energies of the lobes. There is also a significant morphological
difference in that, for the parameters chosen for comparison, the UPP profile
enables more lobe material to escape the centre of the cluster and form two
distinct lobes on either side of the injection cylinder, whereas the
$\beta$-model does not provide sufficient buoyancy to achieve this. In real
clusters, cold material surrounding the cluster core impedes the progression
of jets; in addition, rising bubbles lift material away from the cluster core
which would otherwise have been accreted onto the SMBH and so starve it of
material from which to form the jets. This process would result in the jets
carrying less energy and so reduce the energy carried away by the bubbles.
Such thermostatic control of an AGN is a key element of the feedback loop of
real clusters (McNamara & Nulsen, 2007, 2012; Yang & Reynolds, 2016b), which
regulates the power of the jets; our models have a constant power and so we do
not model this aspect of AGN feedback.

Our simulations have highlighted the role of instabilities in the progression
of the jet. Where the jet is exposed to the ambient medium through the
buoyancy of the lobe having removed it from the core of the AGN, we observe
the growth of KH instabilities; although we believe these have been
exaggerated by our jet injection boundary condition through the formation of a
slower-moving sheath around the jet. While short-wavelength modes of the KH
instability on the surface of the jets are suppressed through resolution
effects, this happens in a similar way for all simulations. The preferential
suppression of the KH instability in the high power jets is hence likely
mainly caused by the shielding by the low-density lobes. This drives the
latter to high aspect ratios, greater than what is generally seen in
observations. Our work highlights the limitations of such numerical
simulations and points to the need to introduce a mechanism to increase the
jet working surface at extended distances, this can be done with the
introduction of a ‘wobble’ such as that used by Perucho et al. (2019, 2022).
Lobes wider than those seen in our simulations can easily be produced using a
precession mechanism, which can be caused by a supermassive secondary or a
misaligned accretion disc (Krause et al., 2019; Horton et al., 2020a; Horton
et al., 2020b, 2023). The latter can, however, also be related to a binary
supermassive black hole (Nealon et al., 2022).

Where our images depart from current observations, this is due to the lack of
non-thermal particles in our model; this means that shock acceleration cannot
be simulated and so the knots and hot spots normally seen are absent. Ground-
breaking work is being carried out by others on this area (see Mukherjee et
al. (2021b) and Yates-Jones et al. (2021)) which have been successful in
reproducing the features of DSA. Similarly, our models do not take into
account radiative losses or spectral ageing, unlike the work of Yates-Jones et
al. (2021).

Future papers will use this UPP atmosphere model to probe the synthetic X-ray
emission from the hot gas at the centre of the cluster and also to further
investigate synthetic radio polarisation images of these clusters, extending
our previous work to these more realistic models.

## Acknowledgements

MJH acknowledges support from the UK STFC [ST/V000624/1]. This work has made
use of the University of Hertfordshire Science and Technology Research
Institute high performance computing facility (https://uhhpc.herts.ac.uk/). We
also thank an anonymous referee for a careful and constructive reading of this
paper which led to significant improvements upon the original text.

## Data Availability

Data available upon request.

## References

  * Antonuccio-Delogu & Silk (2010) Antonuccio-Delogu V., Silk J., 2010, MNRAS, 405, 1303
  * Arnaud et al. (2005) Arnaud M., Pointecouteau E., Pratt G. W., 2005, A&A, 441, 893
  * Arnaud et al. (2010) Arnaud M., Pratt G. W., Piffaretti R., Böhringer H., Croston J. H., Pointecouteau E., 2010, A&A, 517, A92
  * Babyk et al. (2018) Babyk I. V., McNamara B. R., Nulsen P. E. J., Russell H. R., Vantyghem A. N., Hogan M. T., Pulido F. A., 2018, ApJ, 862, 39
  * Barnes et al. (2021) Barnes D. J., Vogelsberger M., Pearce F. A., Pop A.-R., Kannan R., Cao K., Kay S. T., Hernquist L., 2021, MNRAS, 506, 2533
  * Basson & Alexander (2003) Basson J. F., Alexander P., 2003, MNRAS, 339, 353
  * Bell (1978) Bell A. R., 1978, MNRAS, 182, 147
  * Bennett & Sijacki (2022) Bennett J. S., Sijacki D., 2022, MNRAS, 514, 313
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition
  * Blandford & Rees (1974) Blandford R. D., Rees M. J., 1974, MNRAS, 169, 395
  * Blandford et al. (2019) Blandford R., Meier D., Readhead A., 2019, ARA&A, 57, 467
  * Böhringer et al. (2016) Böhringer H., Chon G., Kronberg P. P., 2016, A&A, 596, A22
  * Bonafede et al. (2010) Bonafede A., Feretti L., Murgia M., Govoni F., Giovannini G., Dallacasa D., Dolag K., Taylor G. B., 2010, A&A, 513, A30
  * Bourne & Sijacki (2017) Bourne M. A., Sijacki D., 2017, MNRAS, 472, 4707
  * Bourne & Sijacki (2021) Bourne M. A., Sijacki D., 2021, MNRAS, 506, 488
  * Bourne & Yang (2023) Bourne M. A., Yang H.-Y. K., 2023, Galaxies, 11, 73
  * Bourne et al. (2019) Bourne M. A., Sijacki D., Puchwein E., 2019, MNRAS, 490, 343
  * Carilli & Rawlings (2004) Carilli C. L., Rawlings S., 2004, New A Rev., 48, 979
  * Cavaliere & Fusco-Femiano (1976) Cavaliere A., Fusco-Femiano R., 1976, A&A, 500, 95
  * Cielo et al. (2014) Cielo S., Antonuccio-Delogu V., Macciò A. V., Romeo A. D., Silk J., 2014, MNRAS, 439, 2903
  * Croston et al. (2008) Croston J. H., et al., 2008, A&A, 487, 431
  * Donahue et al. (2006) Donahue M., Horner D. J., Cavagnolo K. W., Voit G. M., 2006, ApJ, 643, 730
  * Eckert et al. (2019) Eckert D., et al., 2019, A&A, 621, A40
  * Ehlert et al. (2018) Ehlert K., Weinberger R., Pfrommer C., Pakmor R., Springel V., 2018, MNRAS, 481, 2878
  * Ehlert et al. (2021) Ehlert K., Weinberger R., Pfrommer C., Springel V., 2021, MNRAS, 503, 1327
  * Ehlert et al. (2023) Ehlert K., Weinberger R., Pfrommer C., Pakmor R., Springel V., 2023, MNRAS, 518, 4622
  * Einstein (1905) Einstein A., 1905, Annalen der Physik, 322, 891
  * English et al. (2016) English W., Hardcastle M. J., Krause M. G. H., 2016, MNRAS, 461, 2025
  * English et al. (2019) English W., Hardcastle M. J., Krause M. G. H., 2019, MNRAS, 490, 5807
  * Fabian (1994) Fabian A. C., 1994, ARA&A, 32, 277
  * Feretti et al. (1995) Feretti L., Dallacasa D., Giovannini G., Tagliani A., 1995, A&A, 302, 680
  * Frenk et al. (1988) Frenk C. S., White S. D. M., Davis M., Efstathiou G., 1988, ApJ, 327, 507
  * Gaibler et al. (2009) Gaibler V., Krause M., Camenzind M., 2009, MNRAS, 400, 1785
  * Gaibler et al. (2011) Gaibler V., Khochfar S., Krause M., 2011, MNRAS, 411, 155
  * Ghirardini et al. (2017) Ghirardini V., Ettori S., Amodeo S., Capasso R., Sereno M., 2017, A&A, 604, A100
  * Ghisellini et al. (2014) Ghisellini G., Tavecchio F., Maraschi L., Celotti A., Sbarrato T., 2014, Nature, 515, 376
  * Gilkis & Soker (2012) Gilkis A., Soker N., 2012, MNRAS, 427, 1482
  * Giri et al. (2022) Giri G., Dubey R. P., Rubinur K., Vaidya B., Kharb P., 2022, MNRAS, 514, 5625
  * Gitti et al. (2012) Gitti M., Brighenti F., McNamara B. R., 2012, Advances in Astronomy, 2012, 950641
  * Godfrey & Shabala (2013) Godfrey L. E. H., Shabala S. S., 2013, ApJ, 767, 12
  * Gopal-Krishna & Wiita (2000) Gopal-Krishna Wiita P. J., 2000, A&A, 363, 507
  * Govoni et al. (2010) Govoni F., et al., 2010, A&A, 522, A105
  * Guo (2015) Guo F., 2015, ApJ, 803, 48
  * Guo (2016) Guo F., 2016, ApJ, 826, 17
  * Gupta et al. (2017) Gupta N., et al., 2017, MNRAS, 467, 3737
  * Hardcastle (2018) Hardcastle M. J., 2018, MNRAS, 475, 2768
  * Hardcastle & Krause (2013) Hardcastle M. J., Krause M. G. H., 2013, MNRAS, 430, 174
  * Hardcastle & Krause (2014) Hardcastle M. J., Krause M. G. H., 2014, MNRAS, 443, 1482
  * Hardcastle & Sakelliou (2004) Hardcastle M. J., Sakelliou I., 2004, MNRAS, 349, 560
  * Hardcastle & Worrall (2000) Hardcastle M. J., Worrall D. M., 2000, MNRAS, 319, 562
  * Hardcastle et al. (1997a) Hardcastle M. J., Alexander P., Pooley G. G., Riley J. M., 1997a, MNRAS, 288, L1
  * Hardcastle et al. (1997b) Hardcastle M. J., Alexander P., Pooley G. G., Riley J. M., 1997b, MNRAS, 288, 859
  * Hardcastle et al. (2005) Hardcastle M. J., Sakelliou I., Worrall D. M., 2005, MNRAS, 359, 1007
  * Hardcastle et al. (2016) Hardcastle M. J., et al., 2016, MNRAS, 455, 3526
  * Harwood et al. (2016) Harwood J. J., et al., 2016, MNRAS, 458, 4443
  * Harwood et al. (2017) Harwood J. J., et al., 2017, MNRAS, 469, 639
  * He et al. (2021) He Y., Mansfield P., Rau M. M., Trac H., Battaglia N., 2021, ApJ, 908, 91
  * Heinz et al. (2006) Heinz S., Brüggen M., Young A., Levesque E., 2006, MNRAS, 373, L65
  * Hillel & Soker (2016) Hillel S., Soker N., 2016, MNRAS, 455, 2139
  * Hillel & Soker (2017) Hillel S., Soker N., 2017, ApJ, 845, 91
  * Hitomi Collaboration et al. (2016) Hitomi Collaboration et al., 2016, Nature, 535, 117
  * Hofmann et al. (2016) Hofmann F., Sanders J. S., Nandra K., Clerc N., Gaspari M., 2016, A&A, 585, A130
  * Horton et al. (2020a) Horton M. A., Hardcastle M. J., Read S. C., Krause M. G. H., 2020a, MNRAS, 493, 3911
  * Horton et al. (2020b) Horton M. A., Krause M. G. H., Hardcastle M. J., 2020b, MNRAS, 499, 5765
  * Horton et al. (2023) Horton M. A., Krause M. G. H., Hardcastle M. J., 2023, MNRAS, 521, 2593
  * Huarte-Espinosa et al. (2011) Huarte-Espinosa M., Krause M., Alexander P., 2011, MNRAS, 418, 1621
  * Jaffe (1980) Jaffe W., 1980, ApJ, 241, 925
  * Jeyakumar et al. (2005) Jeyakumar S., Wiita P. J., Saikia D. J., Hooda J. S., 2005, A&A, 432, 823
  * Kaiser & Alexander (1997) Kaiser C. R., Alexander P., 1997, MNRAS, 286, 215
  * Kardashev (1962) Kardashev N. S., 1962, Soviet Ast., 6, 317
  * King (1962) King I., 1962, AJ, 67, 471
  * Koessl & Mueller (1988) Koessl D., Mueller E., 1988, A&A, 206, 204
  * Komissarov & Porth (2021) Komissarov S., Porth O., 2021, New A Rev., 92, 101610
  * Krause (2005) Krause M., 2005, A&A, 431, 45
  * Krause & Camenzind (2001) Krause M., Camenzind M., 2001, A&A, 380, 789
  * Krause & Camenzind (2003) Krause M., Camenzind M., 2003, New A Rev., 47, 573
  * Krause et al. (2012) Krause M., Alexander P., Riley J., Hopton D., 2012, MNRAS, 427, 3196
  * Krause et al. (2019) Krause M. G. H., et al., 2019, MNRAS, 482, 240
  * Kunz et al. (2011) Kunz M. W., Schekochihin A. A., Cowley S. C., Binney J. J., Sanders J. S., 2011, MNRAS, 410, 2446
  * Laing (1996) Laing R. A., 1996, in Ekers R. D., Fanti C., Padrielli L., eds, Vol. 175, Extragalactic Radio Sources. p. 147
  * Leahy & Williams (1984) Leahy J. P., Williams A. G., 1984, MNRAS, 210, 929
  * Lind et al. (1989) Lind K. R., Payne D. G., Meier D. L., Blandford R. D., 1989, ApJ, 344, 89
  * Longair et al. (1973) Longair M. S., Ryle M., Scheuer P. A. G., 1973, MNRAS, 164, 243
  * Lyskova et al. (2019) Lyskova N., Churazov E., Zhang C., Forman W., Jones C., Dolag K., Roediger E., Sheardown A., 2019, MNRAS, 485, 2922
  * Martí (2019) Martí J.-M., 2019, Galaxies, 7, 24
  * Martizzi et al. (2019) Martizzi D., Quataert E., Faucher-Giguère C.-A., Fielding D., 2019, MNRAS, 483, 2465
  * Massaglia et al. (1992) Massaglia S., Trussoni E., Bodo G., Rossi P., Ferrari A., 1992, A&A, 260, 243
  * Massaglia et al. (1996) Massaglia S., Rossi P., Bodo G., Ferrari A., 1996, Astrophysical Letters and Communications, 34, 295
  * Massaglia et al. (1997) Massaglia S., Bodo G., Rossi P., Ferrari A., 1997, Mem. Soc. Astron. Italiana, 68, 163
  * Massaglia et al. (2016) Massaglia S., Bodo G., Rossi P., Capetti S., Mignone A., 2016, A&A, 596, A12
  * Massaglia et al. (2022) Massaglia S., Bodo G., Rossi P., Capetti A., Mignone A., 2022, A&A, 659, A139
  * Mathews (1971) Mathews W. G., 1971, ApJ, 165, 147
  * Matthews & Scheuer (1990) Matthews A. P., Scheuer P. A. G., 1990, MNRAS, 242, 623
  * McCarthy et al. (1991) McCarthy P. J., van Breugel W., Kapahi V. K., 1991, ApJ, 371, 478
  * McNamara & Nulsen (2007) McNamara B. R., Nulsen P. E. J., 2007, ARA&A, 45, 117
  * McNamara & Nulsen (2012) McNamara B. R., Nulsen P. E. J., 2012, New Journal of Physics, 14, 055023
  * Meece et al. (2017) Meece G. R., Voit G. M., O’Shea B. W., 2017, ApJ, 841, 133
  * Melville et al. (2016) Melville S., Schekochihin A. A., Kunz M. W., 2016, MNRAS, 459, 2701
  * Mendygral et al. (2012) Mendygral P. J., Jones T. W., Dolag K., 2012, ApJ, 750, 166
  * Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228
  * Mignone et al. (2010) Mignone A., Rossi P., Bodo G., Ferrari A., Massaglia S., 2010, MNRAS, 402, 7
  * Miniati (2015) Miniati F., 2015, ApJ, 800, 60
  * Miniati & Beresnyak (2015) Miniati F., Beresnyak A., 2015, Nature, 523, 59
  * Morsony et al. (2010) Morsony B. J., Heinz S., Brüggen M., Ruszkowski M., 2010, MNRAS, 407, 1277
  * Mukherjee et al. (2018) Mukherjee D., Bicknell G. V., Wagner A. Y., Sutherland R. S., Silk J., 2018, MNRAS, 479, 5544
  * Mukherjee et al. (2020) Mukherjee D., Bodo G., Mignone A., Rossi P., Vaidya B., 2020, MNRAS, 499, 681
  * Mukherjee et al. (2021a) Mukherjee D., Bicknell G. V., Wagner A. Y., 2021a, Astronomische Nachrichten, 342, 1140
  * Mukherjee et al. (2021b) Mukherjee D., Bodo G., Rossi P., Mignone A., Vaidya B., 2021b, MNRAS, 505, 2267
  * Mullin et al. (2008) Mullin L. M., Riley J. M., Hardcastle M. J., 2008, MNRAS, 390, 595
  * Nagai et al. (2007) Nagai D., Kravtsov A. V., Vikhlinin A., 2007, ApJ, 668, 1
  * Nandra et al. (2013) Nandra K., et al., 2013, arXiv e-prints, p. arXiv:1306.2307
  * Navarro et al. (1995) Navarro J. F., Frenk C. S., White S. D. M., 1995, MNRAS, 275, 720
  * Nealon et al. (2022) Nealon R., Ragusa E., Gerosa D., Rosotti G., Barbieri R., 2022, MNRAS, 509, 5608
  * Nelson et al. (2014) Nelson K., Lau E. T., Nagai D., Rudd D. H., Yu L., 2014, ApJ, 782, 107
  * Norman et al. (1982) Norman M. L., Winkler K. H. A., Smarr L., Smith M. D., 1982, A&A, 113, 285
  * O’Dea & Saikia (2021) O’Dea C. P., Saikia D. J., 2021, A&A Rev., 29, 3
  * Pakmor & Springel (2013) Pakmor R., Springel V., 2013, MNRAS, 432, 176
  * Pedelty et al. (1989) Pedelty J. A., Rudnick L., McCarthy P. J., Spinrad H., 1989, AJ, 98, 1232
  * Perucho et al. (2011) Perucho M., Quilis V., Martí J.-M., 2011, ApJ, 743, 42
  * Perucho et al. (2014) Perucho M., Martí J.-M., Quilis V., Ricciardelli E., 2014, MNRAS, 445, 1462
  * Perucho et al. (2019) Perucho M., Martí J.-M., Quilis V., 2019, MNRAS, 482, 3718
  * Perucho et al. (2022) Perucho M., Martí J.-M., Quilis V., 2022, MNRAS, 510, 2084
  * Planck Collaboration et al. (2013) Planck Collaboration et al., 2013, A&A, 550, A131
  * Prasad et al. (2015) Prasad D., Sharma P., Babul A., 2015, ApJ, 811, 108
  * Rawlings & Saunders (1991) Rawlings S., Saunders R., 1991, Nature, 349, 138
  * Reynolds et al. (2002) Reynolds C. S., Heinz S., Begelman M. C., 2002, MNRAS, 332, 271
  * Saikia & Salter (1988) Saikia D. J., Salter C. J., 1988, ARA&A, 26, 93
  * Sanders et al. (2011) Sanders J. S., Fabian A. C., Smith R. K., 2011, MNRAS, 410, 1797
  * Scheuer (1974) Scheuer P. A. G., 1974, MNRAS, 166, 513
  * Schuecker et al. (2004) Schuecker P., Finoguenov A., Miniati F., Böhringer H., Briel U. G., 2004, A&A, 426, 387
  * Siegel et al. (2018) Siegel S. R., et al., 2018, ApJ, 861, 71
  * Su et al. (2021) Su K.-Y., et al., 2021, MNRAS, 507, 175
  * Sun et al. (2011) Sun M., Sehgal N., Voit G. M., Donahue M., Jones C., Forman W., Vikhlinin A., Sarazin C., 2011, ApJ, 727, L49
  * Tang & Churazov (2017) Tang X., Churazov E., 2017, MNRAS, 468, 3516
  * Tanner & Weaver (2022) Tanner R., Weaver K. A., 2022, AJ, 163, 134
  * Taub (1948) Taub A. H., 1948, Physical Review, 74, 328
  * Tregillis et al. (2001) Tregillis I. L., Jones T. W., Ryu D., 2001, ApJ, 557, 475
  * Tribble (1991) Tribble P. C., 1991, MNRAS, 253, 147
  * Vazza et al. (2009) Vazza F., Brunetti G., Kritsuk A., Wagner R., Gheller C., Norman M., 2009, A&A, 504, 33
  * Vazza et al. (2012) Vazza F., Roediger E., Brüggen M., 2012, A&A, 544, A103
  * Vazza et al. (2017) Vazza F., Jones T. W., Brüggen M., Brunetti G., Gheller C., Porter D., Ryu D., 2017, MNRAS, 464, 210
  * Vikhlinin et al. (2006) Vikhlinin A., Kravtsov A., Forman W., Jones C., Markevitch M., Murray S. S., Van Speybroeck L., 2006, ApJ, 640, 691
  * Voit et al. (2005) Voit G. M., Kay S. T., Bryan G. L., 2005, MNRAS, 364, 909
  * Wagner & Bicknell (2011) Wagner A. Y., Bicknell G. V., 2011, ApJ, 728, 29
  * Wagner et al. (2012) Wagner A. Y., Bicknell G. V., Umemura M., 2012, ApJ, 757, 136
  * Walker et al. (2015) Walker S. A., Sanders J. S., Fabian A. C., 2015, MNRAS, 453, 3699
  * Weinberger et al. (2017) Weinberger R., Ehlert K., Pfrommer C., Pakmor R., Springel V., 2017, MNRAS, 470, 4530
  * Weinberger et al. (2023) Weinberger R., et al., 2023, MNRAS, 523, 1104
  * Wicker et al. (2022) Wicker R., Douspis M., Salvati L., Aghanim N., 2022, in mm Universe @ NIKA2 - Observing the mm Universe with the NIKA2 Camera. p. 00046 (arXiv:2111.01490), doi:10.1051/epjconf/202225700046
  * Willott et al. (1999) Willott C. J., Rawlings S., Blundell K. M., Lacy M., 1999, MNRAS, 309, 1017
  * Worrall (2009) Worrall D. M., 2009, A&A Rev., 17, 1
  * Yang & Reynolds (2016a) Yang H. Y. K., Reynolds C. S., 2016a, ApJ, 818, 181
  * Yang & Reynolds (2016b) Yang H. Y. K., Reynolds C. S., 2016b, ApJ, 829, 90
  * Yates-Jones et al. (2021) Yates-Jones P. M., Shabala S. S., Krause M. G. H., 2021, MNRAS, 508, 5239
  * Yates-Jones et al. (2022a) Yates-Jones P. M., Shabala S. S., Power C., Krause M. G. H., Hardcastle M. J., Noh Velastín E. A. N. M., Stewart G. S. C., 2022a, arXiv e-prints, p. arXiv:2212.10059
  * Yates-Jones et al. (2022b) Yates-Jones P. M., Turner R. J., Shabala S. S., Krause M. G. H., 2022b, MNRAS, 511, 5225
  * Zanni et al. (2003) Zanni C., Bodo G., Rossi P., Massaglia S., Durbala A., Ferrari A., 2003, A&A, 402, 949

##  Appendix A Cluster magnetic fields

###  A.1 Theory of cluster magnetic fields

The radial scaling of the cluster magnetic field strength with density
distribution was first put forward by Jaffe (1980). Kunz et al. (2011), from
theoretical work on the viscous heating of the ICM suggest a magnetic field
which scales as $B\propto n_{e}^{1/2}$. This result is also recovered by
Bonafede et al. (2010) from rotation measure (RM) observations of the Coma
cluster. This relation is assumed by many authors (e.g. Miniati, 2015; Miniati
& Beresnyak, 2015) and can be expressed as

$\frac{B^{2}}{8\pi}=\eta\frac{3}{2}nk_{B}T$ (15)

where $\eta$ is the energy density ratio between the magnetic field and the
thermal energy density. If the magnetic field strength decreases as the square
root of the thermal electron density, then the gas must remain in overall
hydrostatic equilibrium since the magnetic energy density decreases in
proportion with the gas energy density. Through their work in comparing models
with observations, Böhringer et al. (2016) find that the magnetic energy
density is two orders of magnitude smaller than the thermal energy with $\eta$
(in Equation 15.) falling in the range $5-$7.5\text{\times}{10}^{-3}$$; this
means that the magnetic fields in clusters are considered to be dynamically
unimportant. In this study we will employ a magnetized cluster with a
dynamically unimportant magnetic energy density scaling with the gas energy
density, producing a distribution dependent upon the model chosen for the
density distribution.

###  A.2 Creating the model cluster magnetic field

We define a vector potential of the form $\tilde{A}(k)=A(k)e^{i\theta(k)}$,
where $k$ is the wave-vector ($k^{2}=k_{x}^{2}+k_{y}^{2}+k_{z}^{2}$), $i$ is
the unitary complex number and $A$ and $\theta$ are the vector’s amplitudes
and phases (compare Tribble (1991)). Values for $\theta(k)$ are then drawn
from a uniform random distribution, whereas values for
$\textbf{A}(\textbf{k})$ are drawn from a Rayleigh probability distribution,
given in polar form as:

$P\left(A,\theta\right)\mathrm{d}A\mathrm{d}\theta=\frac{A}{2\pi|A_{k}|^{2}}\exp\left(-\frac{A^{2}}{2|A_{k}|^{2}}\right)\,\mathrm{d}A\mathrm{d}\theta,$
(16)

where $|A_{k}|^{2}\propto k^{-\zeta/2}$. It is then possible to transform
$\tilde{A}(k)$ to real space by employing an inverse Fourier transform. The
three vector potential components are then used to provide the magnetic field
strength (i.e. $B(x)=\nabla\times A(x)$) resulting in an isotropic divergence-
free tangled magnetic field:

$|B_{k}|^{2}=C_{n}^{2}k^{-n}$ (17)

where $n=\zeta-2$ and $C_{n}^{2}$ is the power spectrum normalisation. Huarte-
Espinosa et al. (2011) then employed a Kolmogorov-like 3D turbulent slope
$n=11/3$, which is based upon both theoretical studies of how CMFs evolve
(Vazza et al., 2009) and observations of pressure fluctuations in the Coma
cluster (Schuecker et al., 2004).

Huarte-Espinosa et al. (2011) allowed the magnetized plasma to relax for a
period of time ($118\,$Myr) before injecting their jets to ensure that the
atmosphere was stable. However, here we follow the method introduced in Paper
2 and cut the scale off below $\sim 3$ pixels to avoid injecting high-spatial-
frequency structure into the simulations (this is done by removing the
wavevectors greater than $1/3$ of the Nyquist frequency in Fourier space).
These high frequencies would have been damped out by numerical diffusion
within a short space of time at the start of the run. This, therefore,
provides a more useful baseline at $t=0$ such that jets can be introduced into
a near-stable environment from the very start of the run, maximising valuable
computer-processing time. For consistency we will follow Paper 3 in setting
the field strength to $0.7\,$nT when the the density (in simulation units) is
unity (a realistic value derived from observations of the Coma cluster by
Feretti et al. (1995)). The introduction of the tangled magnetic field in our
model creates a small amount of turbulence in the initially static ICM due to
magnetic tension; the ratio of kinetic energy to thermal energy as a result of
this was measured to be $\sim${10}^{-6}$$ and so our cluster models are, to a
good approximation, relaxed. Our methodology contrasts with the work of some
other researchers who deliberately set out to investigate dynamic cluster
environments (e.g. Ehlert et al., 2018, 2021). We will follow Paper 3 in
increasing the magnetic field energy density in proportion to the gas energy
density (Equation 15) and so $B\propto\sqrt{n_{e}}$. The variation of the
magnetic field strength with distance from the cluster centre can be seen in
Figure 18 where its variable and multi-scaled nature is evident. Plasma
$\beta$ values were calculated for the atmosphere shown in Figure 18 with
typical values falling in the range $${10}^{2}$-${10}^{3}$$; these coincide
with values found in the literature (e.g. Melville et al., 2016; Bourne &
Yang, 2023) and indicate that the magnetic field of our cluster models is not
dynamically important.

Figure 18: Variation of magnetic field strength with distance from the cluster
centre at $t=0$ for a UPP atmosphere with
$M_{500}=3\times${10}^{14}$\,\text{h}^{-1}_{70}M_{\odot}$. Individual colours
are along separate radii from the centre, the black line is the average of
these. All values are normalized to $7\mu$G.

##  Appendix B Simulated Polarimetry

The polarized emission from the lobe can be characterised by means of the
Stokes parameters. Firstly, relativistic aberration needs to be accounted for;
this is where the angle between the light ray and the velocity direction in
the observer’s frame of reference $\theta_{o}$ will be different to that of
the object’s frame $\theta_{s}$ when moving at relativistic speeds. The
formula is (Einstein, 1905)

$\cos\theta_{o}=\frac{\cos\theta_{s}-\beta}{1-\beta\cos\theta_{s}}$ (18)

where $\beta$ is the velocity in units of the speed of light. Therefore, we
define here $B_{x}$ and $B_{y}$ as the components perpendicular to the
aberration-corrected projection axis. The Stokes $I$ (total intensity) and
Stokes $Q$ and $U$ (polarized intensities) parameters are then calculated (in
simulation units) by summing the following relations over the emission volume

$j_{I}=p\left(B_{x}^{2}+B_{y}^{2}\right)^{\frac{1}{2}(\alpha-1)}(B_{x}^{2}+B_{y}^{2})D^{3+\alpha}$
(19) $j_{Q}=\mu
p\left(B_{x}^{2}+B_{y}^{2}\right)^{\frac{1}{2}(\alpha-1)}(B_{x}^{2}-B_{y}^{2})D^{3+\alpha}$
(20) $j_{U}=\mu
p\left(B_{x}^{2}+B_{y}^{2}\right)^{\frac{1}{2}(\alpha-1)}(2B_{x}B_{y})D^{3+\alpha}$
(21)

where $p$ is the local thermal pressure, $\alpha$ is the power-law synchrotron
spectral index, which is taken to be $\alpha=0.5$ corresponding to an electron
energy index $p=2$ and $\mu$ is the maximum fractional polarization for a
given spectral index: for $\alpha=0.5$: $\mu=(\alpha+1)/(\alpha+5/3)=0.69$. D
is the Doppler factor, given by

$D=\frac{1}{\gamma(1-\beta\cos(\theta))},$ (22)

where $\gamma$ is the Lorentz factor and $\theta$ is the angle between the
projection vector and the velocity vector of the cell. The scaling factor to
convert from simulation units to SI units is given as
$j_{0}=$3.133\text{\times}{10}^{31}$\,$WHz-1sr-1, where a constant observing
frequency of $151\,$MHz has been assumed. The Stokes parameters can then be
used to calculate the linear polarisation $\Pi$

$\Pi=\frac{\sqrt{Q^{2}+U^{2}}}{I}.$ (23)

Given the magnetic field is perpendicular to the electric field, then the
observable direction of the magnetic field (projected onto the sky)
$\theta_{B}$ is given by

$\theta_{B}=\frac{1}{2}\arctan\left(\frac{U}{Q}\right)+\frac{\pi}{2}.$ (24)