# A Numerical Study of the Impact of Jet Magnetic Topology on Radio Galaxy
Evolution

Yi-Hao Chen,1 and Sebastian Heinz1, and Eric Hooper1
1Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter
Street, Madison, WI 53706, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The propagation of active galactic nucleus jets depends both on the
environment into which they propagate and on their internal structure. To test
the impact that different magnetic topologies have on the observable
properties of radio galaxies on kpc scales, we conducted a series of magneto-
hydrodynamic simulations of jets injected with different magnetic field
configurations propagating into a gaseous atmosphere modeled on the Perseus
cluster. The simulations show that the structure of the field affects the
collimation and propagation of the jets on cluster scales and thus the
morphology of the radio lobes inflated by the jets, due both to magnetic
collimation and the development of dynamical instabilities in jets with
different magnetic topologies. In all cases, the simulations show a distinct
reversal of the sychrotron spectral age gradient in the radio lobes about a
dynamical time after the jets turn off due to large scale circulation inside
the radio lobe, driven primarily by buoyancy, which could provide a way to
constrain the age of radio sources in cluster environments without the need
for detailed spectral modeling and thus constrain the radio mode feedback
efficiency. We suggest a robust diagnostic to search for such age gradients in
multi-frequency radio data.

###### keywords:

galaxies: active — galaxies: clusters: intracluster medium — galaxies: jets —
magnetohydrodynamics (MHD) — methods: numerical

††pubyear: 2022

##  1 Introduction

The physics of relativistic jets launched by compact objects present an
ongoing modeling challenge. It is generally accepted that jets are launched by
extracting the rotational energy from the black hole (Blandford & Znajek,
1977, hereafter BZ) or the accretion disk (Blandford & Payne, 1982) through
the threading of magnetic fields into the event horizon or the disk. The jet
power of the BZ mechanism depends on the magnetic flux and spin of the black
hole, which are both still poorly constrained. Large-scale poloidal fields,
which could be advected by the accretion flow from larger radii (Beckwith et
al., 2009) or generated through a dynamo around the compact object (Tout &
Pringle, 1996; Bugli et al., 2014), are likely necessary to accelerate the
jets.

Twisting of the magnetic fields by the rotation of the black hole or the
accretion disk naturally leads to the generation of toroidal fields
(Tchekhovskoy, 2015). It is thought that, in the vicinity of the black hole,
the flow could be dominated by Poynting flux as in the original BZ mechanism
and modeled by Li et al. (2006) and Guan et al. (2014), in which the magnetic
energy density is much larger than the kinetic and thermal energy. If the
magnetic energy is transferred to the particles within the jets, the particles
can be accelerated to highly relativistic velocity through a large range of
scales (e.g. Heinz & Begelman, 2000; Lyubarsky, 2009). Collimation, which
could by facilitated by the magnetohydrodynamic (MHD) structure of the jets
(Nakamura et al., 2007; McKinney & Blandford, 2009) and/or by external
pressure such as disk winds (Porth & Fendt, 2010), helps the jets remain beam-
like for up to hundreds of kpc.

The large scale nebulae inflated by propagating jets, generally referred to as
radio lobes or radio galaxies, are commonly observed to have morphologies
described within the framework proposed by Fanaroff & Riley (1974) (see also
Banfield et al., 2015). Nearby, resolved sources are seen with rich structures
inside their radio lobes (Leahy et al., 1996; Miley, 1980). The energetic
feedback by these radio sources likely plays an important role in the
evolution of galaxies and their environments (see e.g. Fabian, 2012; McNamara
& Nulsen, 2012, for reviews).

The BZ mechanism requires strong poloidal fields threading the event horizon
of the black hole to accelerate the jets. Beckwith et al. (2008) tried
different magnetic field structures, including poloidal, toroidal, and
quadrupolar in global general relativistic magneto-hydrodynamic (GRMHD)
simulations, and found that only predominantly poloidal magnetic fields can
generate jets. In these jets, the pitch angle between the poloidal and
toroidal magnetic fields is an essential factor during the propagation of the
jets. Bromberg & Tchekhovskoy (2016) point out that the pitch angle depends on
the external medium that confines the jets, and a denser medium will result in
a smaller pitch angle, i.e. stronger toroidal fields.

The rapid advance of GRMHD simulations has enabled some critical insights into
jet acceleration and jet-disk coupling, revealing that jets are almost
unavoidable consequences of accretion onto rapidly rotating black holes, and
that the net magnetic flux through the disk plays a critical role in
regulating accretion and jet formation. In particular, in the magnetically
arrested disk (MAD) framework, the accretion flow can choke if the magnetic
field exceeds a critical value, leading to conditions conducive to
acceleration of highly relativistic jets. Still, the magnetic topology within
the jet on scales of tens of kiloparsec remains largely unknown.

Observationally, Gabuzda et al. (2015) found gradients of Faraday rotation
measure across the jet structure in observations with the Very Long Baseline
Array (VLBA). These gradients suggest the existence of toroidal components of
the magnetic field at least in the central few 10-pc region. However, Faraday
rotation shows only magnetic fields along the line of sight and can not tell
the fields perpendicular to it. Pasetto et al. (2021) interpret the double-
helical emission morphology of the M87 jet as indication of a helical field
topology, though this is an indirect argument, rather than a direct
measurement of the magnetic field orientation.

From simple flux freezing arguments, it might appear reasonable to assume that
the poloidal component of the field decays away rapidly as the jet expands
laterally once it has accelerated to its terminal velocity. However, jet
stability arguments suggest (e.g. Begelman, 1998; Barniol Duran et al., 2017)
that purely toroidal fields might not be consistent with the observed
structure of jets. In addition, to infer the asymptotic field topology still
requires knowledge of the acceleration scale, the injection field, and the
overall lateral expansion factor, which in most jets is not directly
observable and depends on the cross-field force balance of the jet and thus on
the magnetic field itself.

Thus, it is fair to say that the detailed magnetic topology inside jets is not
fully understood. The work presented here aims to investigate what imprint, if
any, the magnetic topology of jets might have on the large scale propagation
of jets and the overall dynamical properties of radio lobes they inflate.
Given the insurmountable scale discrepancy between jet acceleration and large
scale propagation, we do not attempt to simulate the ab-initio acceleration
and collimation process. Instead, we focus on the interaction of jets with
ambient gas, assuming most of the magnetic energy has been converted to
kinetic energy on scales smaller than the inlets of our injection nozzle,
roughly about 100 pc away from the black hole. Although the energy flux is
dominated by the kinetic energy in our simulations, magnetic fields are still
dynamically important within the jet.

Numerical experiments investigating the propagation of AGN jets have
proliferated in recent years, motivated by the realization that feedback by
jets is crucial for the evolution of galaxies and their environments. For
example, Mendygral et al. (2011); Mendygral et al. (2012) performed 3D non-
relativistic MHD simulations of jets, similar to the ones presented here,
however, with higher plasma $\beta$. Weinberger et al. (2017); Bourne &
Sijacki (2021); Smith & Donohoe (2019, 2021) present long-duration
hydrodynamic simulations of cluster-scale feedback, focusing on the long-term
co-evolution of jets and clusters. In a pioneering set of papers, Horton et
al. (2020); Yates-Jones et al. (2021) investigate the propagation of jets in
clusters using relativistic hydrodynamic simulations, alleviating the
shortcomings of sub-relativstic simulations typically employed in this
context, including the work presented here (see §2.2 and §2.4 for a discussion
on the numerical necessity for this choice). In a complementary long-term
study, Hardcastle & Krause (2013, 2014); English et al. (2016, 2019) present
detailed magneto-hydrodynamic simulations of radio galaxies. All of these
simulations approach the overall problem using different numerical methods and
assumptions.

The work presented in this paper builds on and complements this already rich
canon of published work, focusing on a specific question—the role played by
magnetic topology in the inflation of radio lobes in galaxy clusters—rather
than attempting to perform end-to-end simulations of cosmological structure
formation with feedback. Specifically, we focus on jet propagation in the
well-studied nearby Perseus cluster, with the expectation that findings are
transferable to other cool core clusters at a minimum.

The paper is organized as follows: In §2, we describe the methodology and
techniques used in the simulations. In §3, we present the results of the
simulations, including the dynamics and synchrotron spectral properties, §4
presents a discussion of the results, and §5 summarizes our findings.

##  2 Numerical Techniques

The scale-problem in jet modeling is well known: The physical launching radius
of the jets is a few gravitational radii ($r_{g}=2GM/c^{2}=10^{-4}$ pc for a
$10^{9}M_{\odot}$ black hole). The largest radio lobes can, however, extend to
several Mpc (e.g.,J1420-0545 and 3C 236, Machalski et al., 2008; Tremblay et
al., 2010). It is currently impossible to simulate the acceleration of jets
and the inflation of radio lobes together, which span over 10 orders of
magnitude in spatial scales. Thus, studies of jets can be categorized into
several regimes: (a) jet launching, (b) collimation, and (c) propagation, jet-
ambient interaction. Due to the highly nonlinear physical processes and the
complex and dynamic environments, applications of analytical solutions are
usually limited. Thus, numerical simulations are necessary to study these
systems. Here, we focus on simulations that model regimes (b) and (c).

We use the grid-based adaptive mesh-refinement (AMR) magneto-hydrodynamic
(MHD) code FLASH111http://flash.uchicago.edu/ (Fryxell et al., 2000; Dubey et
al., 2009) to perform the simulations. Full 3D simulations are conducted to
model the non-linear nature of MHD. We include simulations of back-to-back
jets injected into both hemispheres. This allows us to model jet-driven large
scale flows across the mid-plane and to use a comparison of the two sides to
estimate the uncertainty in derived quantities due to the stochasticity of
non-linear fluid mechanics. The simulations under-gird and extend the work
presented in Chen et al. (2019) and Heinrich et al. (2021) and follow the same
overall setup.

We use the unsplit staggered mesh (USM) scheme implemented in FLASH by Lee
(2013) for solving the MHD equations. The USM algorithm in FLASH is a finite-
volume Godunov scheme employing the constrained transport method that
preserves the divergence-free nature of magnetic field to machine precision.
We use the hybrid Riemann solver, which combines the HLL and HLLD solvers for
high accuracy and stability, 3rd- order reconstruction and the mc slope
limiter. This solver is manifestly non-relativistic, the limitations and
implications of which we discuss further below.

The initial conditions are described in detail in §2.1. The simulation box is
$1\times 0.5\times 0.5$ Mpc with $16\times 8\times 8$ base blocks and
hydrostatic diode boundary conditions (hydrostatic-f2+nvdiode). Each AMR block
contains $8^{3}$ cells. The minimum cell size is 30 pc under 12 levels of
adaptive refinement. The jet nozzle diameter is resolved by 2 blocks, or 16
cells. The refinement criteria include the native second derivative based
estimator and also a self-defined momentum-based condition that ensures that
the jets are always covered at the highest level of refinement. We restrict
the maximum refinement level for the ICM further away from the origin such
that the number of blocks is always larger than 32 to the central plane (x-y
plane) and larger than 16 to the the jet axis (z-axis).

The jets are active at a constant jet power of $L_{\rm jet}\approx
10^{45}\,{\rm ergs/s}$ for 10 Myr. The maximum refinement level (used to
resolve the jets) is reduced after the jets are turned off, and the simulation
is allowed to evolve passively for another 91 Myrs to study the long term
evolution of radio lobes. While the choice of duty cycle is arbitrary, the
size scale of the simulated radio lobes corresponds reasonably well with the
observed cavity sizes in the centers of cool core clusters.

Most of the visualization and analysis work is performed with the python-based
software package yt222http://yt-project.org/ (Turk et al., 2011). The 3D
rendering is carried out in VisIt (Childs et al., 2012).

Table 1: List of variables used in the simulation. Variable | Symbol | Value | Unit
---|---|---|---
jet power | $L$ | $1\times 10^{45}$ | erg/s
jet velocity | $v_{j}$ | $3\times 10^{9}$ | cm/s
jet gamma | $\gamma_{j}$ | $4/3$ |
jet density | $\rho_{j}$ | $1.73\times 10^{-26}$ | g/cm3
jet nozzle radius | $r_{\text{jet}}=r_{2}$ | $7.5\times 10^{20}$ | cm
|  | 243 | pc
jet magnetic field | $B$ | $1.7\times 10^{-4}$ | gauss
jet plasma beta | $\beta_{p}$ | 1 |
jet internal Mach number | $M_{j}$ | 10 |
jet toroidal to poloidal field ratio | $h$ |  $\infty$, 1, 0 |
mean molecular weight | $\mu$ | 0.61 |
ICM core density | $\rho_{0}$ | $9.6\times 10^{-26}$ | g/cm3
ICM density profile beta | $\beta$ | 0.53 |
ICM density profile core radius | $r_{c}$ | 26 | kpc
ICM core temperature | $T_{\text{core}}$ | 3.0 | keV
|  | $3.5\times 10^{7}$ | K
ICM outer temperature | $T_{\text{out}}$ | 6.4 | keV
|  | $7.4\times 10^{7}$ | K
ICM temperature profile core radius | $r_{\text{c,T}}$ | 60 | kpc
ICM gamma | $\gamma_{\text{ICM}}$ | $5/3$ |

###  2.1 Initial Conditions

We simulate jet propagation into a non-magnetized cluster initially in
hydrostatic equilibrium and inject highly magnetized
($\beta_{p}=p_{\text{gas}}/p_{B}\sim 1$) jets at the center of the cluster.
Initial density and temperature profiles are set according to the
observational constraints of the Perseus Cluster derived by Zhuravleva et al.
(2015):

The density profile follows a spherically symmetric $\beta$-model:

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

The temperature profile is set such that the core of the cluster is cooler
than the outskirts appropriate to cool core clusters like Perseus:

$T(r)=\frac{T_{\text{out}}[1+(\frac{r}{r_{\text{c,T}}})^{3}]}{\frac{T_{\text{out}}}{T_{\text{core}}}+(\frac{r}{r_{\text{c,T}}})^{3}}.$
(2)

The respective numerical value of each variable can be found in Table 1.

The jet and the ICM fluids follow adiabatic equations of state with
$\gamma_{j}=4/3$ and $\gamma_{\text{ICM}}=5/3$, respectively. The
gravitational potential of the cluster, dominated by dark matter, is held
constant throughout the simulation (thus eliminating the need to solve
Poisson’s equation for self-gravity).

###  2.2 Injection of the Jets

Following Heinz et al. (2006); Morsony et al. (2010), jets are injected into
the grid through a cylindrical “nozzle” at the center of the simulation
domain. The total energy output is governed by the choice of nozzle radius,
velocity $v_{j}$, jet mass density $\rho_{j}$, Mach number $M_{j}$, and
toroidal magnetic field strength $B_{\phi}$ inside the nozzle. The power of
the jets $L$ can be expressed as the sum of the kinetic energy flux, the
enthalpy flux, and the Poynting flux, which can be expressed as the integral
over the faces of the cylindrical nozzle

$L=2\int
v_{j}\left(\frac{1}{2}\rho_{j}v_{j}^{2}+\frac{\gamma_{j}}{\gamma_{j}-1}P_{j}+\frac{B_{\phi}^{2}}{8\pi}\right)\cdot\mathop{}\\!\mathrm{d}{A},$
(3)

where the factor of 2 accounts for two faces of the nozzle. The velocity
profile $v_{\rm z}(r)$ of the jet across the nozzle is constant within the
nozzle radius ($r_{2}$ in Fig. 1) and smoothly goes to zero (at
$r_{\text{out}}$) in a tapering feather region set to be 1/4 of the nozzle
radius for numerical stability:

$v_{j}(r)=\begin{dcases}v_{0}&\text{if }0<r\leq r_{2}\\\
\frac{v_{0}}{2}\left(1+\cos(\pi\frac{r-r_{2}}{r_{\text{out}}-r_{2}})\right)&\text{if
}r_{2}<r\leq r_{\text{out}}\\\ 0&\text{if }r>r_{\text{out}}\end{dcases}$ (4)

We use the jet’s internal sonic Mach number $M_{j}=10$ to parameterize the
ratio between density and pressure given a fixed fluid velocity

$M_{j}=\frac{v_{j}}{c_{s}}=v\sqrt{\frac{\rho_{j}}{\gamma_{j}P_{j}}}.$ (5)

Given that the jet is highly supersonic, the total flux is dominated by the
kinetic energy term, with the other fluxes contributing only at the percent
level.

The plasma $\beta_{p}$, or the ratio of thermal to magnetic pressure

$\beta_{p}=\frac{p_{j}}{p_{B}}$ (6)

defines the strength of the magnetic field inside the jet and is set to 1 for
all magnetized jets. The magnetic pressure within the jet includes both
toroidal ($B_{\phi}$) and poloidal ($B_{z}$) fields inside the nozzle, such
that the magnetic pressure is

$p_{B}=\frac{B^{2}}{8\pi}=\frac{B_{\phi}^{2}+B_{z}^{2}}{8\pi}$ (7)

Note that enforcing $M_{j}=10$ and $\beta_{p}=1$ to be identical for all MHD
simulations—a choice made to facilitate direct comparison—introduces very
small (less than $\sim$1%) variations in total jet power between simulations,
which we consider immaterial.

The nozzle size is then set by the requirement that the jet be in net lateral
pressure balance with its environment at injection. If this were not the case
and the nozzle size were too large, strong re-collimation shocks would form,
crushing the jet and quickly leading to numerical dissipation of some of the
internal toroidal field, thus affecting our ability to control the internal
field structure of the jet (if the jet nozzle were set too small,
computational requirements would unnecessarily increase). The overall maximum
resolution is set by our goal to resolve the jet with at least 16 cells across
both transverse directions of the nozzle.

Together with the condition of $M=10$, the resolution constraint informs our
choice of jet injection velocity, based on computational feasibility: While a
speed as close to $v_{\rm j}=c$ would be ideal (even in non-relativistic
simulations), the computational demand is proportional to the number of
computational cells $N_{3D}$ and the Courant time step $\Delta t_{c}$,
depending on resolution $\Delta R$ as $N_{3D}\propto\Delta R^{-3}$ and $\Delta
t\propto\Delta R/v_{\rm j}$, respectively. The resolution itself depends on
the injection velocity as $\Delta R\propto v_{\rm j}^{-1/2}$, and thus the
total computation demand $N_{\rm compute}\propto N_{3D}\times N_{\rm
step}\propto N_{3D}/\Delta t_{\rm c}$ depends on injection velocity as $N_{\rm
compute}\propto v_{\rm j}^{3}$, in the worst-case scenario. Realistically, the
number of cells in an AMR simulation grows less rapidly than described above,
leading to a somewhat less dramatic dependence of $N_{\rm
compute}\sim\ln{(\Delta R)}\Delta R^{2}/v\propto\ln{v_{\rm j}}\sim\Delta
R^{2}/v_{\rm j}\propto v_{\rm j}^{2}$.

We thus chose $v=0.1c$ for computational feasibility, given the available
resources. We discuss our mitigation strategies for this approximation further
in §2.4. Clearly, future experiments at full relativistic speeds would be
desirable once computationally feasible.

###  2.3 Magnetic Field Injection

We set up 3 different cases of magnetized jets with toroidal, helical, and
poloidal fields in the jets (as depicted in Fig. 2), along with a non-
magnetized case (hereafter referred to as the hydro case) as a control case.

Figure 1: Radial profiles of the magnetic fields (upper and middle panel) and
velocity (lower panel) inside and outside the nozzle. The inner tapering
radius $r_{1}$ is set to be 1/2 of the nozzle radius. The tapering feather
($r_{\text{out}}-r_{2}$, $r_{3}-r_{\text{out}}$, and $r_{\text{pol}-r_{4}}$ )
is set to be 1/4 of the nozzle radius. The poloidal field ($B_{z}$) closes
outside of the nozzle, where we can see negative values. The $B_{z}$ outside
of the nozzle is set to be 10% of the field in the nozzle and $r_{4}$ is the
corresponding radius to close the field loop.

We use a polynomial expression for $B_{\phi}$ to construct a smooth toroidal
magnetic field inside the jet, which is constant ($B_{\phi 0}$) in most parts
of the jet and approaches zero at the axis $r=0$ (to avoid diverging hoop
stresses) and the boundary $r=r_{\text{jet}}=r_{\rm out}$. We split $B_{\phi}$
into 3 regions: 0 to $r_{1}$, $r_{1}$ to $r_{2}$, and $r_{2}$ to
$r_{\text{out}}=r_{\rm jet}$. The polynomials are constructed such that the
derivatives and thus B-fields are continuous at the boundaries. The following
functions are used to set up the toroidal field:

$B_{\phi}=B_{\phi
0}\begin{dcases}-\frac{2r^{3}}{r_{1}^{3}}+\frac{3r^{2}}{r_{1}^{2}}&\text{if
}0<r<r_{1}\\\ 1&\text{if }r_{1}<r<r_{2}\\\
-\frac{2(r_{\text{out}}-r)^{3}}{(r_{\text{out}}-r_{2})^{3}}+\frac{3(r_{\text{out}}-r)^{2}}{(r_{\text{out}}-r_{2})^{2}}&\text{if
}r_{2}<r<r_{\text{out}}\end{dcases}$ (8)

The poloidal field $B_{z}$ is set to be constant ($B_{z0}$) inside the jet and
closed on field lines outside the nozzle such that the magnetic loop is
closed. The initial $B_{z}$ outside of the nozzle is set at 10% of the
strength inside a region that guarantees zero net flux. We construct the
functional form of $B_{z}$ to ensure smooth derivatives at the transitions.
The resulting profile of the magnetic fields is shown in Fig. 1.

Figure 2: Schematic drawings of the injection magnetic topology of the jet:
toroidal (_left_), poloidal (_middle_), and helical (_right_), which is a
weighted sum of the other two cases. The total magnetic field strength inside
the region $r_{1}\leq r\leq r_{2}$ is set to be the same for all 3 cases. The
shaded regions indicate the size of the jet nozzle. Note that the cartoon
shows the poloidal field in the plane of paper, but the actual field is axi-
symmetric.

In practice, we use the magnetic vector potential $\mathbf{A}$ to prescribe
the magnetic field to ensure we maintain a divergence-free field.

$\displaystyle A_{r}=$ $\displaystyle 0$ (9a) $\displaystyle A_{\phi}=$
$\displaystyle\frac{1}{r}\int rB_{z}\mathop{}\\!\mathrm{d}r$ (9b)
$\displaystyle A_{z}=$ $\displaystyle\int B_{\phi}\mathop{}\\!\mathrm{d}r$
(9c)

In order to replenish the toroidal field as the fluid moves out of the nozzle,
we calculate the advection rate and add the appropriate amount of flux to
maintain the same magnetic field strength throughout the active period of the
jet. This field injection region is downstream of the hydrodynamic nozzle
(since the fluid properties inside the nozzle are not updated but instead
serve as an interior inflow boundary condition).

The helicity parameter $h$ is used to control the ratio between toroidal and
poloidal components inside the jet nozzle

$h=\frac{B_{\phi 0}}{B_{z0}}.$ (10)

We set $h=\infty$ for toroidal case, $h=1$ for helical case, and $h=0$ for
poloidal case. The total magnetic field pressure is set to be the same for all
3 cases, such that

$(B_{\phi}^{\text{toroidal}})^{2}=(B_{\phi}^{\text{helical}})^{2}+(B_{z}^{\text{helical}})^{2}=(B_{z}^{\text{poloidal}})^{2}$
(11)

###  2.4 Jet Jitter

We include a modest, random variation in the jet direction around the mean jet
axis, which is oriented along the $z$-axis of the simulation box. This jitter
is imposed for two reasons:

First, it is known from observations that jet axes change direction on
timescales comparable to the jet travel time (e.g. Nulsen et al., 2002; Young
et al., 2005). This may be due to bona fide jet precession, but even in the
absence of an ordered, secular variation in jet axis, dynamical instabilities
in the accretion disks and at the base of the jets are likely to cause
moderate changes in orientation (Pringle, 1996; Heinz & Begelman, 2000).
Extreme re-orientation between episodes of jet activities could even result in
X-shaped radio sources that have been observed (Dennett-Thorpe et al., 2002;
Roberts et al., 2018). This, however, is not what we are modeling in this
work.

Second, analytic and numerical models of uni-directional jets do not generally
reproduce the well-known shape of radio galaxies, leading to much larger
aspect ratios of cocoons and radio lobes than observed (e.g. Norman et al.,
1982; Marti et al., 1994), insufficient to couple the AGN energy to the
central cluster (Vernaleo & Reynolds, 2006). This led Scheuer (1974, 1982) to
propose the so-called dentist drill effect, which spreads jet thrust over a
much larger cone than that prescribed by the jet’s opening angle. Several
numerical studies also show that precession of the jets is necessary to match
the observed morphology of the radio emission (see e.g. Heinz et al., 2006;
Nawaz et al., 2016).

This is particularly important in simulations like ours that are not fully
relativistic: Relativistic jets with a velocity close to $c$ have an energy-
to-momentum ratio of about $c$, while in our simulation the speed of the jets
is set to be 0.1 $c$ and has a correspondingly smaller energy-to-momentum
ratio. The excess momentum leads to an elongated shape that does not fit
observations (see Fig. 3), but can be remedied by distributing this momentum
over a larger solid angle.

Figure 3: Slices showing density through the central plane in simulations with
constant jet direction (non-precessing) at 2.06 and 4.12 Myrs (top and bottom
four panels, respectively). In each age group, from top to bottom: pure
toroidal, helical, poloidal, no magnetic field. One can see the effect of the
magnetic field on the propagation of jets. The toroidal component of the
magnetic field imposes hoop stress on the flow and helps collimate the jets,
while the poloidal component increases the lateral pressure and makes the jet
expand, but stabilizes it against kink instability. See text for more
discussion.

We implement this effect by allowing the direction of the jet to drift in a
random walk prescription, limited within a set opening angle. One can see the
jets are bent in Fig. 4 due to the introduction of this jitter. The resulting
motion of the jet axis throughout the simulation is shown in Fig. 5. The
opening angle of this cone is modest—the average angle between the jet and the
$z$-axis is about 8 degrees, broadly similar to the typical change one might
expect due to dynamical instabilities within the jet (as explored in a rich
body of literature on dynamical stability of relativistic jets; see, e.g.,
Bodo et al. 2019 and references therein) and potential changes in jet
direction from either disk precession or change in the net angular momentum of
accreted material (Bicknell & Begelman, 1996; Caproni & Abraham, 2004; Erlund
et al., 2007; Martí-Vidal et al., 2011; O’Sullivan et al., 2012; Babul et al.,
2013; Britzen et al., 2018; Dominik et al., 2021) and what is required to
explain observed aspect ratios of radio lobes (Scheuer, 1982; Heinz et al.,
2006).

It is important to point out that the jitter pattern is fully randomized
between simulations analyzed in this paper, i.e., each simulation with jitter
follows a statistically similar but different random pattern. We also
performed a set of simulations with identical dither patters (otherwise
following the same range of magnetic topologies) to verify that our results
are robust. We find that material differences arise only during the first 3
Myrs.

To isolate effects unrelated to this jitter in our simulations, we also
performed a control set of simulations of non-precessing jets with otherwise
identical parameters, but run for a shorter time of 5Myrs. Results from these
control runs are presented in §3.1 and §4.1. Unless otherwise mentioned, our
discussion throughout the paper focuses on the set of simulations including
jitter.

###  2.5 Solving the Transport Equation for Tracer Particles

The synchrotron intensity of a synthetic radio map generated from a simulation
frame depends on (a) the number density of electrons (and potentially
positrons), which is derived from the energy density in the simulations,
assuming the internal energy of the non-thermal plasma is dominated by the
relativistic leptons, and (b) the strength and orientation of the magnetic
fields. To obtain the spectra of the synchrotron emission, one further needs
an estimate of (c) the energy distribution of the electrons. There are several
ways to account for (c), including assuming a simple power-law (e.g.
Hardcastle & Krause, 2014; English et al., 2016), grid-based electron
population transport model, for example in Del Zanna et al. (2006); Mendygral
et al. (2012), or Lagrangian particle-based electron transport in (Mimica et
al., 2009; Vaidya et al., 2016, 2018).

Figure 4: Projection of particle $\gamma_{\text{cut}}$, which corresponds to
the cut-off energy of the electron distribution, at 10 Myr. Each particle,
which represents an ensemble of electrons, is accelerated in the termination
shocks, where it starts from a power-law distribution with infinitely large
cut-off energy. As particles propagate with the fluid, their cut-off energy is
changed due to synchrotron radiation and adiabatic expansion/compression.

To this end, we implement a simple yet accurate numerical framework to
calculate synchrotron emission in our simulations, employing the Lagrangian
tracer particle framework in FLASH (Dubey et al., 2012), whereby the radiative
and adiabatic expansion/compression losses/gains of relativistic
electrons/positrons are tracked by individual computational Lagrangian tracer
_particles_ (as opposed to the Eulerian jet tracer fluid that tracks the
fraction of jet fluid in each cell, which our simulations also employ). Note
that the term _particle_ here refers to the numerical perspective of a
Lagrangian fluid element that moves spatially during the simulation. Each
_particle_ represents an ensemble of non-thermal leptons that are modeled in
the spatial proximity of the location of the tracer particle.

Following Bicknell & Melrose (1982); Coleman & Bicknell (1988), we implement
the well-known characteristic solution to the particle transport equations
subject to radiative losses and adiabatic compression.

We assume the synchrotron radiation from the radio lobes is optically thin to
self-absorption, and the radiation energy is negligible compared to the
thermal energy of the gas. This is consistent with most extended radio sources
observed, which have steep spectral indices (see e.g. Laing et al., 1983).

Particles are injected at the cross-sections of the jet nozzle and are assumed
to initially represent a power-law distribution of relativistic electrons

$f(\gamma,t=t_{0})=N_{0}\gamma^{-p}$ (12)

for $\gamma\geq\gamma_{\rm min}=10$. However, any injection spectrum can in
principle be modeled this way in post-processing.

While $p$ is a free parameter, we choose $p=2$ throughout, given the optically
thin synchrotron spectral index of
$\alpha=-d\log{F_{\nu}}/d\log{nu}=(p-1)/2\sim 0.5$ for typical young radio
sources. The normalization $N_{0}$ is then related to the non-thermal particle
pressure $P$, as

$N_{0}=\frac{3P}{m_{\rm e}c^{3}\ln{\gamma_{\rm max}/\gamma_{\rm min}}}$ (13)

Given the very weak logarithmic dependence on the limits of the normalization
on $\gamma_{\rm min}$ and $\gamma_{\rm max}$, we set $\gamma_{\rm min}=10$,
while we track the evolution of $\gamma_{\rm max}$ as described below.

For each _particle_ , the mapping of the initial distribution to the time-
evolved distribution is accomplished using the characteristic solution for
particles injected with energy $\gamma_{0}$ as a function of time:

$\gamma_{\text{c}}(\gamma_{0},t)=\frac{\left(\frac{n(t)}{n_{0}}\right)^{\frac{1}{3}}}{1/\gamma_{0}+\int_{t_{0}}^{t}A(t^{\prime})\,(\frac{n(t^{\prime})}{n_{0}})^{\frac{1}{3}}\mathop{}\\!\mathrm{d}t^{\prime}}$
(14)

where $n(t)$ is the number density of electrons at time $t$, $n_{0}$ is the
initial number density of _leptons_ , and
$A=\frac{4\sigma_{T}}{3m_{e}c}\left(U_{B}+U_{\rm CMB}\right)$ is the radiative
energy loss rate, neglecting synchrotron-self-Compton losses (which are
insignificant for electrons in diffuse synchrotron emitting regions like radio
lobes). For an injection spectrum of $p=2$, the solution reduces to a simple
power-law with the same index but a sharp cutoff at

$\gamma_{\rm
cut}=\lim_{\gamma_{0}\rightarrow\infty}\gamma_{c}(\gamma_{0},t)=\frac{\left(\frac{n(t)}{n_{0}}\right)^{\frac{1}{3}}}{\int_{t_{0}}^{t}A(t^{\prime})\,(\frac{n(t^{\prime})}{n_{0}})^{\frac{1}{3}}\mathop{}\\!\mathrm{d}t^{\prime}}$
(15)

In the simulations, we have to choose a time $t_{0}$ to begin the integration,
i.e., a reference time when the electrons are accelerated to have a power-law
energy distribution. The acceleration mechanism of particles inside jets is
still not fully understood; the site of the acceleration could be located
inside the jet (Bicknell & Melrose, 1982) and/or at terminal hotspots, and the
mechanism could be, e.g., magnetic reconnection (Romanova & Lovelace, 1992;
Sironi & Spitkovsky, 2014; Sironi et al., 2015; Guo et al., 2015) or diffusive
shock acceleration (Bell, 1978; Blandford & Ostriker, 1978; Drury, 1983;
Caprioli & Spitkovsky, 2014; Park et al., 2015).

We choose to begin the cooling calculation when the fluid element the particle
travels in has dissipated most of its kinetic energy, identified as the time
at which the bulk velocity of the _particles_ first drops below half of the
jet injection velocity.

This approach requires significantly less computational effort and memory than
solving the full transport equations numerically (Mendygral et al., 2012;
Vaidya et al., 2018). More importantly, the existence of a closed, analytic
solution to the transport equation makes this approach more accurate than a
finite-difference solution of the same equation.

Figure 5: The pattern of the jet pointing direction as implemented in the
simulation. The pointing undergoes a random walk process but is constrained by
a maximal angle that can be set in the simulation. Red crosses label the
position and time of the pointing since the launch of the jet. The labels on
the radial-axis are in units of degrees, from 0 to 20∘.

###  2.6 Synchrotron Emissivity

We map the synchrotron spectrum of individual tracer particles onto the grid
and produce synthesized radio observations of the simulated data in post-
processing by projecting along a given line-of-sight. The synchrotron spectra
of the grid cells are determined following the standard expressions for the
integrated synchrotron emission coefficient (e.g. Rybicki & Lightman, 1979)
for an isotropic pitch angle distribution.

Given the sharp cutoff to the particle distribution the shape of the
synchrotron spectrum for each tracer particle is well-described by a simple
power-law with an exponential cutoff, such that the emission coefficient in
each cell is

$j_{\nu}\approx\frac{3\left(B\sin{\alpha}\right)^{3/2}e^{7/2}N_{0}}{8\sqrt{\pi}c^{5/2}m_{\rm
e}^{3/2}}\nu^{-1/2}C_{0}e^{-\left(\nu/\nu_{\rm max}\right)}$ (16)

with $\alpha$ being the angle between the magnetic field vector and the line
of sight such that $B_{\perp}=B\sin{\alpha}$ is the in-sky component of the
B-field, and, for $p=2$,

$C_{0}=\frac{4\sqrt{2}}{3}\Gamma\left(2\frac{1}{12}\right)\Gamma\left(\frac{5}{12}\right)\approx
4.16$ (17)

The cutoff frequency is

$\nu_{\rm max}=\frac{3\gamma_{\rm max}^{2}eB\sin{\alpha}}{4\pi m_{\rm e}c}$
(18)

in which $\gamma_{\text{max}}$ is evaluated for each tracer particle using eq.
(15); the normalization $N_{0}$ is calculated using eq. (13); and pressure $P$
and magnetic field $B$ are taken from the simulation grid.

In practice, we calculate the value of $j_{\nu}$ in each cell through a
weighted average of the nearest 10 particles, weighted by inverse square
distance from the cell center. This is necessary because not every
computational cell contains at least one tracer particle.

The emission coefficient can then be integrated along the line-of-sight to
calculate the synthesized synchrotron intensity maps, which are presented in
Section 3.3. Given the fact that the emission clock of each tracer particle is
started after the jet fluid has dissipated its kinetic energy, we remove any
emission from the jet fluid itself (i.e., removing any computational cell with
velocity in excess of half the jet speed) and focus our analysis of the
synchrotron emission on the radio lobes.

##  3 Results

A motivating goal of our simulations is to investigate how different topology
of the magnetic fields in the jets affect the jet-ambient medium interaction,
including the collimation, propagation, and the imprints on the radio lobes,
and to probe which physical effect is responsible for any observed
differences. We thus compare the outcome of simulations with different
injected magnetic topology to each other and to a purely hydrodynamic control
simulation.

###  3.1 Non-Precessing Jets

We perform a set of control simulations without explicit precession to serve
as our calibration and baseline cases. These non-precessing simulations have
the same configuration as our precessing simulations, but without the jitter
described in §2.4. The simulations were run for a shorter duration of 5 Myrs.

Snapshots of density slices are shown in Fig. 3. We can see that each jet
propagates at different speed into the external medium. Although the jet
velocity is the same across all simulations, the head of the jet is influenced
by the internal dynamics of magnetic fields and thus has different ability to
propagate.

We determine the propagation velocity of the jet head in the simulations by
measuring the endpoint of the jet tracer fluid in the direction of the jet
axis (z-axis). This is achieved by creating profiles of the jet fluid mass
fraction against the z-coordinate. We then determine the location of the jet
head to be where this fraction $f_{\rm W}$ exceeds $10^{-6}$. We analyze both
sides (+z and -z) and the difference between them can be considered as the
numerical uncertainty for the analysis.

Fig. 6 shows the results of this analysis. The edges of the shaded regions
represent the location of the jet heads. We find that the toroidal jet
propagates the fastest at the beginning (before 1 Myr) followed by the helical
and hydro jets, while the poloidal jet is the slowest. Later, the hydro jet
becomes slightly faster than all others between 1 and 3 Myr, but slows down
after about 3 Myr. Finally, after 4 Myr, the poloidal jet remains the slowest,
while the hydro jet becomes slower than the helical and toroidal jets. The
hydro jet propagates at a velocity comparable to the poloidal jet.

Figure 6: Propagation distances of the jet heads for the 4 non-precessing
jets. The edges of the filled region represent the longest extension of the
jet fluid in the +z and -z directions for simulations without precession.
Dotted and dashed lines: linear propagation rate for comparison. Note that the
jet fluid is set to a velocity of 0.1 $c$ $\sim$ 30 kpc/Myr in our
simulations.

Analytically, the propagation velocity of the jet head can be estimated by
balancing the ram pressure in the frame of the jet head once the density ratio
and the bulk velocity in the jets is known (see e.g. Norman & Winkler, 1985).
The velocity of the jet head can be written as

$v_{h}=\frac{\sqrt{\eta}}{1+\sqrt{\eta}}v_{j},$ (19)

where $v_{j}$ is the jet fluid velocity and $\eta=\rho_{j}/\rho_{e}$ is the
ratio of the jet fluid density and environmental density. In our cases, the
jet velocity is fixed to be 0.1 $c$ ($\sim$ 30 kpc/Myr), and the environmental
density profile is the same for all simulations. The variation of the
propagation velocity thus arises from the density of the jet. Although the
density of the nozzle is the same for all simulations, the different
configurations of magnetic fields cause the jet to expand or self-collimate at
different rates. Fig. 7 shows the density profile along the jet axis at about
4 Myr. Two reference lines of constant head velocity $v_{h}$ are plotted in
Fig. 6: 8 and 10 kpc/Myr for density ratio $\eta$ of 0.12 and 0.23
respectively.

Even though in all simulations the jets begin with the same density in the
nozzle, the magnetic structure determines the later evolution of the jet. The
figure shows that the jet in the toroidal case (green line) maintains the
highest mean density as it propagates, while the jet in the poloidal case
drops to very low density as soon as the fluid leaves the nozzle.

Figure 7: Averaged density profile within 0.25 kpc from the axis of the non-
precessing jets at 4 Myr. All simulations begin with the same density in the
nozzle, but the different configurations of the magnetic fields give rise to
different density profiles after the jet fluid leaves the nozzle. This causes
the different propagation velocity of the jet heads in our four cases.

Based purely on the lateral force balance of a steady jet, one would thus
expect that the toroidal jet should propagate fastest due to the hoop stresses
imposed by the toroidal magnetic fields on the flow, which collimate the jet,
compress the flow, and increase the density, leading to a higher momentum
flux, thus enabling the fluid to penetrate faster in the ICM. While this is
true initially, we find that the toroidal jet slows down to speeds below that
of the helical case.

On the other hand, the poloidal magnetic fields act as additional lateral
pressure that expands the jet and reduces the density, making it slower to
propagate, which suggests that the poloidal case should have the slowest
propagation velocity.

However, the high density and degree of collimation of the toroidal case do
not consistently make the toroidal jet propagate the fastest. As seen in Fig.
6 and the right panel of Fig. 3, the helical jet propagates as fast as the
toroidal jet at a later time ($\gtrsim$ 3 Myr), despite the lower toroidal
field and the added lateral pressure due to the poloidal field. The reasons
for this are described in detail in Section 4.1.

###  3.2 Jet Propagation with Jitter

When jitter is implemented (Section 2.4), the shapes of the cavities look very
different from the non-precessing jets. The cavities created by the jets
resemble the lobes of typical FR-II radio galaxies (Fanaroff & Riley, 1974, ;
some typical examples include Cygnus A, Centaurus A, Fornax A and Hydra A),
whereas in the case without jitter, the cavities/lobes are very skinny, unlike
those of actual radio galaxies. In Fig. 8, we show the central slices of the
four simulations at about 10 Myr. In all cases, the jets have been
continuously active for 10 Myr, after which the jet is turned off and the
lobes are allowed to evolve passively for another 91 Myrs.

The lobe shapes of the toroidal and helical cases are noticeably different
from the poloidal and hydro cases. When toroidal fields are included in the
jet (both toroidal and helical cases), the cavities are more elongated, with
larger aspect ratios, compared to those without toroidal fields (poloidal and
hydro cases), with narrower waists in the equatorial region of the radio
galaxy.

Figure 8: 2

9Slices showing the density through the central plane from four simulations
with different magnetic field topology in the jets at 10 Myr. Top left: pure
toroidal, top right: helical, bottom-left: poloidal, bottom-right: no magnetic
field.

To quantify the differences in propagation, we once again measure the extent
of the lobes from our simulations, similar to the non-precessing case. The
location of the jet heads is plotted in Fig. 9. At the very early stage ($<$ 1
Myr), we see the same sequence of toroidal $>$ hydro $\gtrsim$ helical $>$
poloidal in lobe sizes. However, after 1 Myr the toroidal case slows down to
speeds comparable to the helical case, and the hydro case slows down to speeds
comparable to the poloidal case. At 10 Myr, the time shown in Fig. 8, they
form two distinct groups with toroidal $\sim$ helical $>$ hydro $\sim$
toroidal. As pointed out in §2.4, the simulations with four different magnetic
topologies follow different randomized jitter patterns. In order to verify
that the differences in lobe propagation are unaffected by this randomization,
we repeated this analysis using re-simulations with identical dither pattern.
The only material difference between these two sets of simulations appear
during the first 3 Myrs. Because our analysis of the simulations with jitter
focuses on the later stages of the simulations, we conclude that this
difference at early times does not affect our conclusions in any meaningful
way.

Figure 9: Evolution of the lobe extent in the four precessing simulations.
Edges of the shaded area represent the furthest extents of the jet fluid in +z
and -z directions, representing two sides (upper and lower) of the jets and
can be interpreted as numerical uncertainties in the simulations. The jets are
turned off at 10 Myr. Note that both axes are in log-scale.

The behavior can be understood in the context of the non-precessing control
runs in Section 3.1. The toroidal and helical cases, given the action of
magnetic hoop stress, can collimate better and carry more ram pressure during
the propagation, resulting in the cavities being more elongated. On the other
hand, the poloidal and hydro cases lack any internal collimation mechanism and
propagate more slowly, after an initial period during which lateral pressure
balance is established.

However, the pure toroidal case once again does not propagate faster than the
helical case. This is in contradiction to the understanding that the
propagation velocity is determined by the collimation of the jet and thus by
the strength of toroidal fields. We have seen similar results in the non-
precessing case, in which the propagation of the pure toroidal case slows down
to be comparable to the helical case.

Similarly, the hydrodynamic case approaches the propagation speed of the
poloidal field case, despite the absence of the additional lateral pressure
resulting from the poloidal field. We attribute these deviations from the
naive expectations to the development of dynamical instabilities, which we
will discuss in §4.1.

###  3.3 Synchrotron Properties

With the tracer particle implementation outlined above, we can synthesize the
synchrotron radiation at various frequencies of different evolutionary stages
for different magnetic injection topologies. We show a sequencce of
synchrotron emission maps at 150 MHz in Fig. 10. In the early stages, due to
the precession of the jets, the shape of the radio image is spiky with various
episodes of the jet rapidly propagating into the ambient medium as jitter
changes the jet direction. The emission is brightest at the endpoints of the
spikes, where the jet fluid runs into the denser medium and decelerates. These
locations correspond roughly to the hotspots seen in many of the FRII objects
(note again that emission from the jet itself is not included in these maps).

At later time, after 10 Myr when the jets are turned off, the shape of the
radio lobes becomes roughly spherical, resembling bubbles/cavities seen in
cool core clusters. As the plasma ages, the synchrotron intensity drops due to
radiative cooling and adiabatic expansion.

At even later times, these inflated bubbles rise further and form vortex
rings, as has been discussed extensively in the literature on radio sources in
cool core clusters (see e.g. Churazov et al., 2001; Gardini, 2007; Werner et
al., 2010; Churazov et al., 2013). If viewed from an inclined line-of-sight
(Fig. 11), a synchrotron emitting torus is visible at 100 Myr in each
simulation. As proposed by Churazov et al. (2001), the torus is formed due to
the low-density hot gas rising in the stratified cluster atmosphere, subject
to Rayleigh-Taylor instability.

The velocity field in such a vortex ring reflects large scale, organized
rotation about the axis of the torus. The shape of the radio torus resembles,
for example, the morphology of Fornax A (Fomalont et al., 1989; Anderson et
al., 2018), which has passively evolved for at least 100 Myr and does not show
any hotspots, and is also reminiscent of the ”radio ear” seen in the 90cm
image of Virgo A (Owen et al., 2000). This structure is common among all three
field-congurations, though less pronounced in the poloidal case.

Figure 10: Evolution of synchrotron emission at different times in the
simulations as seen at 90 degree inclination of viewing angle to jet axis.
_Left_ column: toroidal magnetic fields; _middle_ column; helical magnetic
fields; _right_ column: poloidal magnetic fields. Figure 11: Same as Fig. 10
but viewed from a 30-degree inclination angle relative to the jet axis. Notice
that the radio emitting plasma has a torus shape in all 3 cases that are
obvious in an inclined angle.

To derive a measure of spectral aging, we calculate the projected intensity at
two frequencies based on the _particle_ aging model described in Section 2.5
and 2.6. The spectral index between two frequencies is then
$\alpha=\log(I_{\nu_{1}}/I_{\nu_{2}})/\log(\nu_{1}/\nu_{2})$. The results are
shown in Figs. 12 and 13. A gradient from the flatter spectra at the heads to
the steeper spectra at the tails is apparent in the figure when the jet is
still active ($<$ 10Myr).

The synchrotron cooling tracer particle implementation adopted for these
figures assumes that the electron distribution is a fresh power-law at the
working surface/termination shock, where the bulk velocity of the plasma drops
and the bulk of the jet kinetic energy is dissipated in a strong shock,
motivating our assumption that the synchrotron emitting plasma is accelerated
in these termination shocks. This results in a young and relatively flatter
spectrum ($\alpha\sim 0.5$) at the hotspots. Thus, when the jets are still
active, the tips of the lobes generally show flatter spectra compared to the
lower regions closer to the AGN sources. Note again that we explicitly filter
out the emission from the jets. In actual observations, the removal of the
active jet component can be crucial to revealing this gradient, especially in
FR-I sources, where the inner lobe emission is dominated by the jet.

After leaving the hotspots, the plasma starts to cool, creating a steeper
spectrum toward the inner parts of the lobes.

Interestingly, after the jets turn off, as the radio lobe evolves, we find
that the initial outside-in age gradient is not just diluted by the internal
motion of the lobe plasma, but in fact reversed. This is due to the vortex
formed during the rising of the bubble. The rotation of the vortex transports
the young plasma down to the outside face of the vortex, while old plasma
rises faster in the center due to the overall process of vortex formation.
This is most apparent in the lowest two rows in Figs. 12 and 13. This reversal
of the synchrotron age gradient could be useful when searching for and
analyzing fossil radio plasma in X-ray cavities of galaxy clusters. We will
discuss this in more detail in §4.3.

The radio lobes in the poloidal case, compared to helical and toroidal cases,
have relatively flatter spectral indices among the 3 magnetized cases. This
can be seen in Fig. 12, where the poloidal case has generally redder false
color at every stage. This effect is most likely caused by the overall weaker
magnetic field in the poloidal field case, which can be seen in Fig. 15. The
cooling rates in the poloidal case are thus relatively slower, yielding
flatter spectra when compared to the other 2 cases at the same time. On the
other hand, the helical case has the steepest spectral index because of the
stronger magnetic fields. Note that although we inject more magnetic energy in
the toroidal case, it does not have the strongest magnetic fields because of
the onset of the kink instability and the resulting numerical reconnection
discussed in Section 4.1.

Figure 12: Evolution of synchrotron spectral index at different times in the
simulations. The spectral index $\alpha$ is calculated as
$\alpha=\log(I_{\nu_{1}}/I_{\nu_{2}})/\log(\nu_{1}/\nu_{2})$. In this figure,
$\nu_{1}$ = 1400 MHz and $\nu_{2}$ = 150 MHz. _Left_ column: toroidal magnetic
fields; _middle_ column; helical magnetic fields; _right_ column: poloidal
column. Figure 13: Same as Fig. 12 but viewed from 30-degree inclination
relative to the jet axis

##  4 Discussion

###  4.1 Magnetic Field Strength and the Development of Kink Instability
within the Jet

We find that the helical jet with jitter overtakes the purely toroidal
geometry to become the fastest propagating jet head at later times. As we show
below, the helical case also maintains the highest total magnetic energy among
all three magnetized cases, while the poloidal case has the lowest. Finally,
we saw that the hydrodynamic case slows down to speeds comparable to the
poloidal case despite the lack of the additional lateral magnetic pressure. We
understand these outcomes in the context of different dynamical instabilities
the simulated jets are subject to, and the resulting differences both in jet
propagation (which is a physical effect) and numerical dissipation (which is,
at root, an un-physical effect that is unavoidable in any grid-based
simulation, but sometimes similar to turbulent dissipation encountered in
astrophysical systems).

Magnetic energy injection occurs differently for the different cases: For the
poloidal case, the initial dipole field is stretched and amplified and no
additional magnetic flux is injected in the simulation domain. On the other
hand, we inject additional toroidal magnetic flux into the domain for the
toroidal and helical cases that is constantly transported out of the nozzle.

Due to the lateral expansion of the jets, the toroidal ($\phi$) component
scales as $1/r$ while the poloidal ($z$) component scales as $1/r^{2}$. For
initially comparable fields in toroidal and poloidal components, once the jet
expands, the toroidal component will dominate the overall magnetic field in
the jet.

Figure 14: 3D rendering of a simulation snapshot with helical magnetic fields
in the jets. At left, slices show the temperature map (hot in blue and cold in
red). At right, streamlines of magnetic fields represent the magnetic
structure inside the lobe. Although the magnetic fields are highly ordered in
the jet, they become mostly tangled in the lobe.

In the lobes, magnetic fields become mostly tangled due to the turbulent
motions. This can be seen in the rendered image in Fig. 14 as well as the
evolution of different magnetic components in Fig. 15. In order to investigate
the evolution of the magnetic field, we measure the magnetic energy of the
$(r,z,\phi)$ components in cylindrical coordinates, with the $z$-axis parallel
to the mean jet axis. In the non-precessing cases (denoted “np“), the $\phi$
component is maintained at a high level for both toroidal and helical cases,
while in the precessing cases, the $\phi$ component is no longer dominant
energetically. Rather, the different components are comparable to each other,
especially around and after 10 Myr when the jet is turned off as shown in the
bottom panel of Fig. 15. This is likely due to the turbulent nature of the
flow within the radio lobes, which efficiently mixes the different modes.

Figure 15: Evolution of the magnetic energy and different components of the
magnetic energy in the simulations. Different colors denote different magnetic
configurations of the injection nozzle of each simulation. Solid lines
indicate the total magnetic energy, while the filled charts show the
percentage of the $r$, $z$, and $\phi$ components. _Top panel_ : non-
precessing simulations (denoted “np“). Although in toroidal case more magnetic
flux is injected in the domain, at later times the amount of magnetic energy
is surpassed by the helical case. The poloidal case remains the lowest in
magnetic energy. _Bottom panel_ : precessing jets active for 10 Myr. The
magnetic energy decreases after the jet is off due to (mostly numerical)
reconnection in the lobes. Note that at later times (after 10Myr), the $z$
component increases, which is likely caused by the buoyantly rising bubble
stretching the magnetic fields vertically.

The long-term behavior of the helical and the purely toroidal configuration
require further discussion. We find that the counter-intuitive behavior
described above can be explained by the action of the kink instability, which
can also account for the field decrease in the toroidal case. This can be
clearly seen in Fig. 16. The toroidal ($\phi$) field in the toroidal case,
although initially stronger than in the helical case, suffers from dynamical
disruption, effectively widening the jet cross section, and subsequent
numerical dissipation of the magnetic field loops due to the kink modes that
become non-linear along the jet.

Appl et al. (2000) pointed out that current-driven instabilities are important
in astrophysical jets and carried out linear stability analysis of the growth
rates. They found that the growth rate depends on the radius of the jet, the
Alfvén velocity, as well as the magnetic pitch $P=rB_{z}/B_{\phi}$. When
normalized by the Alfvén crossing time ($r_{\text{jet}}/v_{A}$), the growth
rate is inversely proportional to the pitch or third power of the pitch
depending on small or large pitch regimes. Larger pitch leads to a lower
growth rate, i.e. the kink instability is suppressed. Although in our
simulations the pitch profile is different from the profiles analyzed in Appl
et al. (2000), we can still compare our instability growth rate to their
calculations.

Figure 16: Central slices of the magnetic fields in the x-direction
(perpendicular to the image plane). In the toroidal case, kink instabilities
start to develop beyond 15 and 20 kpc at 2.06 and 4.12 Myr (top and bottom
three panels, respectively). Due to the limited resolution in our simulations,
the instabilities result in magnetic energy being converted to thermal energy.
Thus, it reduces the jet head propagation velocity of the toroidal case at
later times.

First, we consider the toroidal case, which should have a magnetic pitch close
to zero and develop instabilities very quickly. However, we do not see the
growth of the instability until the jet reaches beyond 18 kpc (Fig. 16). This
can be accounted for by the numerical resolution required to resolve the
unstable modes in our simulations. Roughly speaking, the minimum wavelength
that can be resolved and thus grow in the simulation is about 4 cell-widths,
which is about a quarter of the jet diameter. The maximum wavenumber that can
be resolved is then $k_{\text{max}}=2\pi/\lambda_{\text{min}}\sim
12/r_{\text{jet}}$. The growth time of this unstable wavenumber is about
$3r_{\text{jet}}/v_{A}$ (c.f. Fig. 3 in Appl et al. (2000)) or $\sim 0.2$ Myr
given $r_{\text{jet}}\sim 240$ pc and $v_{A}=0.012c$ in our simulations. The
jet fluid is moving at 0.1 $c$ in our simulations, which gives an instability
e-folding length scale of 6 kpc. The instability becomes visible at about 3
times this length scale.

Second, we consider the helical case, in which the jet has magnetic pitch of
order of 1 ($P/r_{\text{jet}}\simeq 1$). From Fig. 3 in Appl et al. (2000), we
can infer that the maximum growth rate is about 0.1 ($\Gamma
r_{\text{jet}}/v_{A}\sim 0.1$) and thus the e-folding growth time is of order
10 Alfvén crossing times ($t=1/\Gamma=10\,r_{\text{jet}}/v_{A}$), which is
roughly 0.8 Myr in our cases. Considering the jet fluid is 0.1 $c$, this
timescale corresponds to a length scale of 24 kpc. We would expect the
instability to become apparent at 3 or 4 times of this length scale. This can
explain why there is no visible instability for the helical jet in Fig. 16.

This is similar to the findings of previous works (Nakamura et al., 2007; Guan
et al., 2014; Bromberg & Tchekhovskoy, 2016; Barniol Duran et al., 2017) that
if the ratio between toroidal and poloidal fields is large, internal kink
modes will develop and disrupt the jets. Due to the limit of our numerical
resolution, the kink instability drives numerical reconnection at the core of
the jet, which suppresses the magnetic field strength in the toroidal case
below that of the helical case.

In Bromberg & Tchekhovskoy (2016); Barniol Duran et al. (2017), the authors
attribute the development of the kink instability to the external medium
density profile since the pitch angle of their jets depends strongly on the
external medium. They argue that the jets propagating in denser external
medium (headed jets) suffer from the kink instability, while the jets
propagating in the less dense medium or a previously evacuated funnel
(headless jets) can maintain their stability. Various factors, e.g. the spin
of the central black hole, can change the magnetic pitch. In this work, we
confirm that even in the same density profile, different magnetic pitches can
still result in varying growth rates of the kink instability.

Finally, in the hydrodynamic case, in which magnetic fields are absent, we see
the development of Kelvin-Helmholtz (KH) modes at the boundaries of the
cavity. The magnetized jets and lobes suppress the KH instability at the
shearing boundaries. This is consistent with previous findings that magnetic
fields can suppress the instabilities (see e.g. Ryu et al., 2000; Jones & De
Young, 2005).

###  4.2 Comparison to Cygnus A

Given the obvious visual similarities, it is instructive to compare our
synthesized radio maps with Cygnus A, the brightest and nearest powerful FR II
radio source (Carilli & Barthel, 1996). Although our simulations are not
specifically modeled to resemble the Cygnus A system, the resulting radio maps
have a lot of similarities to the observed properties of the proto-typical
FRII object, including the morphology as well as spectral index distribution
(McKean et al., 2016).

In Fig. 10, single or multiple hotspots are visible when the jet is active
(before 10 Myr). The misalignments between jet-axis and hotspots are natural
results of the jitter implementation and mimic the two hotspots seen in Cygnus
A (Pyrzas et al., 2015). While active jets are seen in Cyg A, we compare only
the lobes with our synthesized image since we leave out most of the jet
emission. At 50 Myr, more filamentary structures are developed, similar to
Cygnus A. We note that the estimated age of the Cygnus A ranges from few Myr
(Carilli et al., 1991) to 100 Myr due to the large uncertainties in the
magnetic fields. Differences in structure could also arise from the different
duty cycles of the source.

Our spectral index maps at 10 Myr show a strong gradient from the hotspot to
the tail of the lobe. The head of the lobe, where the plasma is freshly
injected, has a flat spectrum with spectral index $\alpha$ close to 0.5, while
the spectral index steepens toward the tail with $\alpha\sim 2$, where most of
the plasma is old. This is consistent with the detailed observations of Cygnus
A (Carilli et al., 1991).

Clearly, the morphology of the jets simulated here differs from the typical
FRI-type morphology of cool-core cluster radio sources like Perseus A and
Virgo A. It worth emphasizing that, while the true nature of the difference
between FRI and FRII sources is still unsettled, it is likely rooted in
physics we did not include in our simulations (e.g., proper special
relativistic MHD with bulk Lorentz factors of order 5 or higher, mass loading
by entrainment, e.g., Bicknell & Begelman 1996). However, the goal of our
simulations is not a full reconstruction of a particular system, but rather a
broader investigation of the impact of different field topologies on the
properties of radio galaxies.

###  4.3 Spectral index gradient evolution

In Section 3.3, we see a clear spectral index gradient while the jets are
still active when the sources are observed close to perpendicular to the jet
axis. During the active jet phase (shown at the time of 9.5 Myr in the
relevant figures), we see a flatter spectral index/younger plasma located at
the tips of the lobe, while the older plasma is at the tails close to the
central plane. This gradient of the spectral index has been seen in many
observations of FR-II sources (e.g. Steenbrugge et al., 2010; Shulevski et
al., 2017; Savini et al., 2018). While it is more difficult to establish such
global trends in FR-I sources, given the subtraction of the jet emission from
the diffuse emission from the inner radio lobe, it has been observed in some
FR-I sources (e.g. Kolokythas et al., 2015).

At a later time after the jet turns off (roughly at 50 Myr as shown in the
figures), we find an intriguing reversal of the age gradient in the spectral
index map. The relatively younger plasma is now at smaller radii of the lobe,
while the outermost parts of the lobes appear to be older, as shown in Fig.
12.

If the viewing direction is closer to the jet axis, the spectral index is
flatter in the center and gets steeper closer to the edges of the lobe while
the jet is active, as shown in the third row (9.5 Myr) of Fig. 13. However,
once the jets shut off, the distribution of the spectral index becomes again
reversed – steeper at the center and relatively flatter closer to the edges,
as in the 4th row (56.7 Myr) of Fig. 13.

This reversal can be understood as a consequence of rotation within the rising
bubble. Our radio sources resemble typical FR-II sources, in which the radio
emission is brighter away from the central engine. The fresh, younger plasma
appears in the vicinity of the hotspots which are located at the farthest ends
of the radio lobes. While it is still not clear what acceleration mechanism
powers FR-II sources — whether the relativistic electrons are accelerated in
the terminal shocks, i.e. the hotspots, or in situ in the jets —, those
sources produce a spectral index (or age) gradient from flatter (younger) at
the farther ends of the lobes to steeper (older) closer to the central source.

Two effects contribute to this rotation:

  * •

First, the injection of the jet creates a backflow toward the tail of the
lobe, which has been pointed out and extensively discussed by numerous authors
(e.g., Norman et al., 1982; Antonuccio-Delogu & Silk, 2010; Cielo et al.,
2017). Depending on the size and location of the cavity, backflows might reach
back to the equatorial plane where the supermassive black hole is located, or
circulate within the cavity. In our simulations, the helical and toroidal
cases have cavities away from the central plane, while the poloidal and hydro
cases show more puffy cavities extending to the equatorial plane (Fig. 8), due
to the differences in internal collimation between these cases. However, in
all cases, a backflow exists, driven by the pressure gradient within the
cocoon. Figure 17 shows the net circulation generated within the radio lobes
of the helical simulation 1 Myr after the jet terminates, at 11Myrs.

  * •

Second, the buoyantly rising bubbles naturally evolve into vortex rings that
drive the global rotation of the plasma in the bubble about the ring (e.g.
Churazov et al., 2001). The sense of rotation is such that the plasma moves
away from the central source (down the gravitational potential) along the axis
of the vortex (i.e., the mean jet direction), and downwards along the outside
of the vortex ring. That is, the rotation twists the plasma in a double-right-
handed way. Because of the ongoing buoyant rise of the lobes/bubbles, the
vortex z-velocity still has positive net values in most of the lobe, however,
in the frame of the bubble, the gas is moving downward. Figure 18 shows the
importance of buoyancy forces compared through what we term the force balance
$FB$, defined as

$FB\equiv\frac{\left|\nabla P\right|-\left|\rho\vec{g}\right|}{\left|\nabla
P\right|+\left|\rho\vec{g}\right|}$ (20)

Positive and negative values of $FB$ indicate regions where pressure forces or
gravity (buoyancy) forces dominate the momentum equation, respectively.

Figure 17: Large scale circulation within the radio lobes of the helical field
simulation. Shown in each panel: Average z-velocity, weighted by the jet
tracer fluid mass fraction for 11, 21, 51, and 101Myrs from top to bottom,
each showing upwelling along the central spine of the lobe and down-flow along
the sides. Initially, this backflow motion is driven by the pressure-gradient
within the cocoon. At later times, the circulation is driven by buyancy. Both
effects contribute to the reversal of the age gradient seen at later times in
the simulation. Note that the net z-velocity of the lobe is positive;
velocities that are downward in the frame of the lobe still have positive net
z-velocity in the frame of the observer. The limits on the color bar were
chosen arbitrarily to optimize the dynamic range of each figure.

The result of this circulation within the bubble is a reversal of the age
gradient that can be seen at roughly one vortex-turnover time after the jet
shuts of, as shown in Fig. 19.

This suggests an interesting observational diagnostic, as detecting a reversed
age gradient could be used as a dynamical clock. In our simulations, the
vortex-turnover time can be estimated using the buoyancy or sound crossing
time of the lobe cavity,

$t_{\rm vortex}\sim t_{\rm sound}\sim\frac{2R_{\rm lobe}}{c_{\rm s}}\sim
35\,{\rm Myrs}\frac{R_{\rm lobe}}{15\,{\rm kpc}}\frac{900\,{\rm km/s}}{c_{\rm
s}}$ (21)

which indicates that the vortex should have formed at around the 45Myr mark in
our simulations (roughly 35 Myrs after jet turn-off). Because this is a purely
dynamical effect, it should be robust against uncertainties in the
normalization of the equipartition fraction $\eta$ of the magnetic field (both
in the initial conditions and the effects of numerical dissipation), though
the optimal spectral region sensitive to measure age gradients will depend on
$\eta$.

However, it is important to point out that this effect might still operate in
the absence of any buoyancy, at least in part, given the internal circulation
driven within the lobe. In particular, Ruszkowski et al. (2007) point out that
external magnetic fields in the ICM can stabilize buoyant bubbles against
disruption, possibly delaying the formation of a buoyant vortex ring. Because
both buoyancy and cocoon dynamics act in concert, we expect the reversal of
this age gradient to be robust, however, re-simulations will be required that
include a magnetized ICM and probe different levels of stratification of the
environment to verify this statement.

Figure 18: Slice of the net buoyancy force balance corresponding to the
velocity panels in Fig. 17. The figure shows that the initial circulation is
not driven by buoyancy, as the entire cocoon is strongly overpressured and
expanding. However, by 50Myrs, the lobes are in approximate pressure
equilibrium with the ICM on average and buoyancy is the driving mechanism for
large scale motion.

Confirming the predicted reversal may be observationally challenging.
Currently, most observed active FR-II sources exhibit age gradients from the
head to the tail since most of the detected sources are higher luminosity
objects whose powering jets are still active or have been active recently (and
have thus not had time to develop buoyancy-driven vortex rings). Once the jet
stops, the radio luminosity quickly drops at least an order of magnitude
within a few Myr, which makes detecting the predicted reversal of the age
gradient more challenging (since this is a subtle effect, as can be seen in
Fig. 12. Although these low-surface brightness radio sources are brighter at
low frequencies, the inherently lower resolution at lower frequency might make
the detection of the reversing effect more difficult.

To derive a more robust diagnostic, we measure the cumulative flux fraction
along the jet

$CF_{\nu}(z)=\frac{\int_{-z_{\rm min}}^{z}\int_{y_{\rm min}}^{y_{\rm
max}}I_{y,z^{\prime}}dydz^{\prime}}{\int_{-z_{\rm min}}^{z_{\rm
max}}\int_{y_{\rm min}}^{y_{\rm max}}I_{y,z}dydz}$ (22)

as a function of distance along the jet at 150MHz and 1.4GHz (such that each
is a function between 0 and 1).

Given the observed age gradient, we would expect these $CF_{\nu}(z)$ curves to
show a difference in shape depending on whether the emission is more or less
centrally concentrated. In Fig. 20, we then plot the difference $\Delta
CF\equiv CF_{150MHz}-CF_{1.4GHz}$ between these cumulative flux fractions as a
function of $CF_{150MHz}$ for the three different simulations and at four
different times. At early times (12.5Myrs), all three simulations show an
S-shaped curve that indicates an inward age gradient (i.e., more of the high-
frequency emission is concentrated at larger distances from the center), thus,
the curve shows a valley and then a peak from left to right. However, at later
times (50Myrs, 100Myrs), the figures show an overall reversal of this S-shape,
indicating the reversal of the age gradient, with younger/high frequency
emission concentrated closer to the center. While the detailed shapes of these
curves do show some differences, this reversal appears robust and suggests an
observational diagnostic that is both simple to derive from multi-frequency
data and robust against model assumptions.

Because $CF$ is easy to calculate and presents a twice-integrated measure of
intensity, it is relatively robust and substantially less noisy than, for
example, spatial averages of the spectral index as a function of $z$. We
propose using this diagnostic as a measure of spatially resolved spectral
aging.

More work beyond the scope of this publication will be required to establish
the presence of this reversal age gradient and to calibrate $CF$ as a true
universal diagnostic. Still, the promise to quantify the time since injection
of particles by active jets ceased (i.e., the time since the jets shut off)
would be useful in, for example, determining the fraction of time radio
sources inflated in galaxy clusters are actively powered vs. passively
evolving, this providing a new constraint on the duty cycle of the jet.
Conversely, since the synchrotron cooling age depends on $\eta$, establishing
the actual dynamical age of the plasma (one turnover time since the jet turned
off) could be used to constrain $\eta$.

Figure 19: The cutoff frequency of the _particles_ at 3 different times of the
helical case. The size of the points represents the local gas density–bigger
points indicate lower density. While the jet is still active, we can see a
clear gradient from the tip to the tail of the lobe. While the jet is off,
part of the older plasma that was initially at the tail is lifted by the
vortex to the top. This make a _reversed_ age gradient at $\sim$ 60 Myr in our
case.  Figure 20: Difference in cumulative emission along the jet axis between
150MHz and 1.4GHz, plotted against cumulative flux at 150MHz. All three
simulations show a characteristic reversal of the s-shaped curve indicating
the reversal of the age gradient in the radio lobes between 25Myrs and 50Myrs.
The red and green lines in the left-hand panels show a cartoon of the expected
signal of the observed age inversion to highlight the change in curvature. The
features at the center of the emission (around $CF\sim 0.5$) in the poloidal
case are due to the fact that this case maintains relatively strong emission
in the center of the cluster due to the poloidal flux that is preserved even
after the jet turns off.

###  4.4 Caveats and Limitations of Our Simulation

While numerical simulations of jet propagation and radio galaxy evolution are
proliferating in the literature, each work must necessarily make simplifying
assumptions and suffers from unavoidable numerical limitations. It is
therefore important to discuss the limitations and caveats of our work in this
context.

####  4.4.1 Strength of the magnetic fields in the jet and lobes, numerical
dissipation

In our simulations, we do not set the strength of the magnetic fields
explicitly. Instead, we choose the plasma beta $\beta_{p}$ to be of order
unity inside the jet nozzle, which implies that the magnetic pressure is
comparable to the thermal pressure of the gas and is dynamically important
within the jet. The field strength is then determined by our choices of jet
power, internal Mach number, and nozzle cross-section. We note the resulting
magnetic field strength at the jet base ($\sim$ kpc from the black hole) is
about 170 $\mu$G. While it is currently impossible to precisely infer the
strength of the magnetic fields in the jets from observations, O’Sullivan &
Gabuzda (2009) find the magnetic fields in AGN cores are in the range of
$100\sim 300$ $\mu$G at 1 pc and scales like $r^{-1}$. Our setting of magnetic
fields is consistent with their findings.

However, we inject our jets ballistically, i.e., assuming most of the magnetic
and thermal energy has been converted to kinetic energy on scales smaller than
the boundaries of our injection nozzle roughly about 100 pc from the black
hole. This is one possible assumption about the nature of jet launching and
propagation. For example, Sikora et al. (2005) argue that, from observational
data, the conversion of Poynting flux to kinetic energy happens mostly within
$10^{3}$ gravitational radii at parsec scale, which is consistent with our
assumptions. However, other plausible models exist which maintain mostly
poynting-flux dominated jets even on kpc scales (e.g., Colgate et al., 2015),
which are clearly in conflict with the assumptions underlying our simulations.

The relatively small fraction of thermal and magnetic energy in the jets
accounts for the relatively minor differences between different magnetic
topologies in the overall shapes of the sources, since the kinetic power is
approximately two orders of magnitude larger than the thermal and electro-
magnetic power. This also implies that, while the jets are injected with a
plasma beta of order unity, the shock-heating of jet fluid substantially
increases the plasma beta within the lobes, i.e., the field strength in the
radio lobes is well below equipartition.

In our current setup, the surrounding medium is non-magnetized. All the
magnetic flux comes from the injection of the jet. Although this setup is
obviously unphysical, we use our simulations to create a baseline for future
studies in which more complicated factors might affect the dynamics of the
jets and the bubbles. Some research hints that the evolution of bubbles might
be affected more by the surrounding magnetic field than by the field inside
the bubbles (Ruszkowski et al., 2007; Dursi & Pfrommer, 2008). We will leave
the investigation of jets interacting with magnetized surroundings to future
work.

Finally, it is important to address the question of numerical diffusivity,
which is an unavoidable consequence of any Eulerian scheme to integrate the
MHD equations (note once again that FLASH preserves the solenoidal condition
to machine accuracy). The effect of this can be seen in Fig. 15, which shows
that the magnetic energy of the simulations is substantially reduced at late
times.

Some reduction in the total magnetic energy is expected because of the PdV
work done by the magnetic field (as part of its adiabatic evolution). The
reduction in the poloidal case is of order 50%, which is consistent with the
expectation of simple adiabatic behavior, given that the radio lobes expand
roughly by a factor of 2 in radius over this timescale.

However, the reduction by up to 80% is substantially in excess of that
expectation in the helical and toroidal cases. We attribute this reduction to
the numerical reconnection of magnetic energy on small simulation scales in
the radio lobes, which is more prevalent in the toroidal and helical cases,
given the much shorter field-reversal scales. Comparing to the poloidal case,
we conclude that of order 50% of the field energy in the helical and toroidal
simulations is lost to numerical reconnection, beyond the level of turbulent
reconnection we would naturally expect to occur (e.g. Lazarian et al., 2020).
Reducing the magnitude of this effect would be desirable, but would require a
substantially higher resolution, beyond the scope of resources available at
the time of the simulations.

As such, aspects of our simulations subject to the total magnitude of the
magnetic field strength (such as the absolute value of the spectral
break/cutoff frequency, see below, and the absolute value of the synchrotron
intensity) should be considered with these caveats in mind, while qualitative
results, such as the overall spatial distribution of age gradients and the
differences in emission morphology between simulations will be more robust
against uncertainties introduced by numerical dissipation. Because FLASH is an
energy conserving code, any effects of numerical dissipation simply increase
the internal energy of the gas.

####  4.4.2 On Spectral Aging Models

Conventionally, spectral aging models assume a constant magnetic field, which
is usually derived from minimum energy arguments, to calculate the age of the
radio source. However, as seen in our simulations, the magnetic fields that a
particular plasma ensemble experiences along its trajectory are time variable.
In the hotspots where the fresh plasma begins to cool, the magnetic fields are
strongest due to compression at the terminal shocks. When the plasma leaves
the hotspots, the magnetic field strength significantly drops because of the
adiabatic expansion of the plasma as well as reconnection of the highly
tangled field. We note that this latter effect is enhanced in any numerical
simulation because of limited resolution and numerical diffusivity.

As a consequence, the cooling rates in the lobes are much lower compared to
the cooling rates in the hotspots. The evolution of the spectral index depends
on the strength of the magnetic field, due to the reduction in overall field
strength from lobe expansion and numerical reconnection.

The Lagrangian particles begin cooling after they leave the jet. Our current
model does not capture any re-acceleration outside the hotspots. A more
complex framework (e.g. Vaidya et al., 2018) would be necessary to study
various acceleration mechanisms and the resulting synchrotron emission.

##  5 Conclusions

In this paper we present a series of MHD simulations of AGN jets with
different injecting field topology. We set up pure toroidal, pure poloidal,
and helical magnetic fields inside the jets (see Fig. 2). The jets are then
evolved according to the ideal MHD formalism. To this end, we implement an
implicit scheme to integrate the particle transport equation that allows
spectral age modeling of radio plasmas.

Our simulations reproduce the well-known result that non-relativistic, heavy
jets cannot reproduce realistic radio lobe shapes unless a random jitter—i.e.,
a dentist drill effect— is imposed on the orientation of the jet axis.

We find that the presence of toroidal magnetic fields reduces the lateral
expansion of the jet and enables faster propagation of the jet heads. On the
other hand, the presence of poloidal fields inhibits the development of kink
instability. Thus, we find that purely toroidal jets do not propagate faster
than helical jets, despite the stronger effects of self-collimation at equal
field strength. Instead, their propagation speed is comparable to the helical
case, as the kink-instability effectively spreads the momentum flux of the jet
over a larger solid angle than might be expected from a fully self-collimated
jet.

We find that, due to the more efficient collimation, toroidal and helical jets
in our simulations generate more elongated cocoons with more pronounced
hourglass morphology and that, as a result, cavities inflated by these jets
are further away from the central engine compared to those inflated by the
poloidal and hydrodynamic jets in our simulations.

While the initial injection of magnetic energy into the lobes reflects the
topology of the magnetic field in the jet, turbulence within the lobes quickly
re-arranges the field to rough equipartition between the different field
orientations.

We find that, after the jets turn off, the radio lobes in our simulations
develop the well-known rising vortex ring structure of buoyant bubbles in
stratified atmospheres in all cases. When viewed from lines of sight close to
the jet axis, the toroidal structure of the emission is readily observable.

We show that the large scale rotation of the vortex rings leads to the a
reversal of the spectral index gradient after the jet shuts off, due to the
reciprocal motions of the young and old plasma. This suggests that the
observation of spectral age gradients can be used as a dynamical clock. While
this effect is likely most pronounced in FR-II radio galaxies (most similar to
the radio sources produced by our simulations), some FR-I radio galaxies
should show similar reversals. We propose a simple, robust diagnostic to
search for such age gradient signatures.

##  6 Data Availability

The numerical data and data products presented in this paper are hosted on
storage maintained by the Department of Astronomy at the University of
Wisconsin-Madison. Data and products will be shared on reasonable request to
the corresponding author, Sebastian Heinz, at<EMAIL_ADDRESS>
## References

  * Anderson et al. (2018) Anderson C. S., Gaensler B. M., Heald G. H., O’Sullivan S. P., Kaczmarek J. F., Feain I. J., 2018, ApJ, 855, 41
  * Antonuccio-Delogu & Silk (2010) Antonuccio-Delogu V., Silk J., 2010, MNRAS, 405, 1303
  * Appl et al. (2000) Appl S., Lery T., Baty H., 2000, A&A, 355, 818
  * Babul et al. (2013) Babul A., Sharma P., Reynolds C. S., 2013, ApJ, 768, 11
  * Banfield et al. (2015) Banfield J. K., et al., 2015, MNRAS, 453, 2326
  * Barniol Duran et al. (2017) Barniol Duran R., Tchekhovskoy A., Giannios D., 2017, MNRAS, 469, 4957
  * Beckwith et al. (2008) Beckwith K., Hawley J. F., Krolik J. H., 2008, ApJ, 678, 1180
  * Beckwith et al. (2009) Beckwith K., Hawley J. F., Krolik J. H., 2009, ApJ, 707, 428
  * Begelman (1998) Begelman M. C., 1998, ApJ, 493, 291
  * Bell (1978) Bell A. R., 1978, MNRAS, 182, 443
  * Bicknell & Begelman (1996) Bicknell G. V., Begelman M. C., 1996, ApJ, 467, 597
  * Bicknell & Melrose (1982) Bicknell G. V., Melrose D. B., 1982, ApJ, 262, 511
  * Blandford & Ostriker (1978) Blandford R. D., Ostriker J. P., 1978, ApJ, 221, L29
  * Blandford & Payne (1982) Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
  * Bodo et al. (2019) Bodo G., Mamatsashvili G., Rossi P., Mignone A., 2019, MNRAS, 485, 2909
  * Bourne & Sijacki (2021) Bourne M. A., Sijacki D., 2021, MNRAS, 506, 488
  * Britzen et al. (2018) Britzen S., et al., 2018, Monthly Notices of the Royal Astronomical Society, 478, 3199
  * Bromberg & Tchekhovskoy (2016) Bromberg O., Tchekhovskoy A., 2016, MNRAS, 456, 1739
  * Bugli et al. (2014) Bugli M., Del Zanna L., Bucciantini N., 2014, MNRAS, 440, L41
  * Caprioli & Spitkovsky (2014) Caprioli D., Spitkovsky A., 2014, ApJ, 783, 91
  * Caproni & Abraham (2004) Caproni A., Abraham Z., 2004, Monthly Notices of the Royal Astronomical Society, 349, 1218
  * Carilli & Barthel (1996) Carilli C. L., Barthel P. D., 1996, A&ARv, 7, 1
  * Carilli et al. (1991) Carilli C. L., Perley R. A., Dreher J. W., Leahy J. P., 1991, ApJ, 383, 554
  * Chen et al. (2019) Chen Y.-H., Heinz S., Enßlin T. A., 2019, MNRAS, 489, 1939
  * Childs et al. (2012) Childs H., Brugger E., Whitlock B., 2012, in Bethel E., Childs H., Hansen C., eds, ”High Performance Visualization: Enabling Extreme-Scale Scientific Insight.”. Chapman and Hall/CRC, doi:10.1201/b12985
  * Churazov et al. (2001) Churazov E., Brüggen M., Kaiser C. R., Böhringer H., Forman W., 2001, ApJ, 554, 261
  * Churazov et al. (2013) Churazov E., Ruszkowski M., Schekochihin A., 2013, MNRAS, 436, 526
  * Cielo et al. (2017) Cielo S., Antonuccio-Delogu V., Silk J., Romeo A. D., 2017, MNRAS, 467, 4526
  * Coleman & Bicknell (1988) Coleman C. S., Bicknell G. V., 1988, MNRAS, 230, 497
  * Colgate et al. (2015) Colgate S. A., Fowler T. K., Li H., Hooper E. B., McClenaghan J., Lin Z., 2015, The Astrophysical Journal, 813, 136
  * Del Zanna et al. (2006) Del Zanna L., Volpi D., Amato E., Bucciantini N., 2006, A&A, 453, 621
  * Dennett-Thorpe et al. (2002) Dennett-Thorpe J., Scheuer P. A. G., Laing R. A., Bridle A. H., Pooley G. G., Reich W., 2002, MNRAS, 330, 609
  * Dominik et al. (2021) Dominik R. M., Linhoff L., Elsässer D., Rhode W., 2021, Monthly Notices of the Royal Astronomical Society, 503, 5448
  * Drury (1983) Drury L. O., 1983, Reports on Progress in Physics, 46, 973
  * Dubey et al. (2009) Dubey A., Reid L. B., Weide K., Antypas K., Ganapathy M. K., Riley K., Sheeler D., Siegal A., 2009, arXiv e-prints, p. arXiv:0903.4875
  * Dubey et al. (2012) Dubey A., Daley C., ZuHone J., Ricker P. M., Weide K., Graziani C., 2012, ApJS, 201, 27
  * Dursi & Pfrommer (2008) Dursi L. J., Pfrommer C., 2008, ApJ, 677, 993
  * 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
  * Erlund et al. (2007) Erlund M. C., Fabian A. C., Blundell K. M., Moss C., Ballantyne D. R., 2007, MNRAS, 379, 498
  * Fabian (2012) Fabian A. C., 2012, ARA&A, 50, 455
  * Fanaroff & Riley (1974) Fanaroff B. L., Riley J. M., 1974, MNRAS, 167, 31P
  * Fomalont et al. (1989) Fomalont E. B., Ebneter K. A., van Breugel W. J. M., Ekers R. D., 1989, ApJ, 346, L17
  * Fryxell et al. (2000) Fryxell B., et al., 2000, ApJS, 131, 273
  * Gabuzda et al. (2015) Gabuzda D. C., Knuettel S., Reardon B., 2015, MNRAS, 450, 2441
  * Gardini (2007) Gardini A., 2007, A&A, 464, 143
  * Guan et al. (2014) Guan X., Li H., Li S., 2014, ApJ, 781, 48
  * Guo et al. (2015) Guo F., Liu Y.-H., Daughton W., Li H., 2015, ApJ, 806, 167
  * 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
  * Heinrich et al. (2021) Heinrich A., Chen Y., Heinz S., Zhuravleva I., Churazov W., 2021, MNRAS, submitted
  * Heinz & Begelman (2000) Heinz S., Begelman M. C., 2000, ApJ, 535, 104
  * Heinz et al. (2006) Heinz S., Brüggen M., Young A., Levesque E., 2006, MNRAS, 373, L65
  * Horton et al. (2020) Horton M. A., Krause M. G. H., Hardcastle M. J., 2020, MNRAS, 499, 5765
  * Jones & De Young (2005) Jones T. W., De Young D. S., 2005, ApJ, 624, 586
  * Kolokythas et al. (2015) Kolokythas K., O’Sullivan E., Giacintucci S., Raychaudhury S., Ishwara-Chandra C. H., Worrall D. M., Birkinshaw M., 2015, MNRAS, 450, 1732
  * Laing et al. (1983) Laing R. A., Riley J. M., Longair M. S., 1983, MNRAS, 204, 151
  * Lazarian et al. (2020) Lazarian A., Eyink G. L., Jafari A., Kowal G., Li H., Xu S., Vishniac E. T., 2020, Physics of Plasmas, 27, 012305
  * Leahy et al. (1996) Leahy J. P., Bridle A. H., Strom R. G., 1996, in Ekers R. D., Fanti C., Padrielli L., eds, Vol. 175, Extragalactic Radio Sources. p. 157
  * Lee (2013) Lee D., 2013, Journal of Computational Physics, 243, 269
  * Li et al. (2006) Li H., Lapenta G., Finn J. M., Li S., Colgate S. A., 2006, ApJ, 643, 92
  * Lyubarsky (2009) Lyubarsky Y., 2009, ApJ, 698, 1570
  * Machalski et al. (2008) Machalski J., Kozieł-Wierzbowska D., Jamrozy M., Saikia D. J., 2008, ApJ, 679, 149
  * Martí-Vidal et al. (2011) Martí-Vidal I., Marcaide J. M., Alberdi A., Pérez-Torres M. A., Ros E., Guirado J. C., 2011, A&A, 533, A111
  * Marti et al. (1994) Marti J. M., Mueller E., Ibanez J. M., 1994, A&A, 281, L9
  * McKean et al. (2016) McKean J. P., et al., 2016, MNRAS, 463, 3143
  * McKinney & Blandford (2009) McKinney J. C., Blandford R. D., 2009, MNRAS, 394, L126
  * McNamara & Nulsen (2012) McNamara B. R., Nulsen P. E. J., 2012, New Journal of Physics, 14, 055023
  * Mendygral et al. (2011) Mendygral P. J., O’Neill S. M., Jones T. W., 2011, ApJ, 730, 100
  * Mendygral et al. (2012) Mendygral P. J., Jones T. W., Dolag K., 2012, ApJ, 750, 166
  * Miley (1980) Miley G., 1980, ARA&A, 18, 165
  * Mimica et al. (2009) Mimica P., Aloy M. A., Agudo I., Martí J. M., Gómez J. L., Miralles J. A., 2009, ApJ, 696, 1142
  * Morsony et al. (2010) Morsony B. J., Heinz S., Brüggen M., Ruszkowski M., 2010, MNRAS, 407, 1277
  * Nakamura et al. (2007) Nakamura M., Li H., Li S., 2007, ApJ, 656, 721
  * Nawaz et al. (2016) Nawaz M. A., Bicknell G. V., Wagner A. Y., Sutherland R. S., McNamara B. R., 2016, MNRAS, 458, 802
  * Norman & Winkler (1985) Norman M. L., Winkler K.-H. A., 1985, Supersonic jets
  * Norman et al. (1982) Norman M. L., Winkler K. H. A., Smarr L., Smith M. D., 1982, A&A, 113, 285
  * Nulsen et al. (2002) Nulsen P. E. J., David L. P., McNamara B. R., Jones C., Forman W. R., Wise M., 2002, ApJ, 568, 163
  * O’Sullivan & Gabuzda (2009) O’Sullivan S. P., Gabuzda D. C., 2009, MNRAS, 400, 26
  * Owen et al. (2000) Owen F. N., Eilek J. A., Kassim N. E., 2000, ApJ, 543, 611
  * O’Sullivan et al. (2012) O’Sullivan E., et al., 2012, Monthly Notices of the Royal Astronomical Society, 424, 2971
  * Park et al. (2015) Park J., Caprioli D., Spitkovsky A., 2015, Phys. Rev. Lett., 114, 085003
  * Pasetto et al. (2021) Pasetto A., et al., 2021, ApJ, 923, L5
  * Porth & Fendt (2010) Porth O., Fendt C., 2010, ApJ, 709, 1100
  * Pringle (1996) Pringle J. E., 1996, MNRAS, 281, 357
  * Pyrzas et al. (2015) Pyrzas S., Steenbrugge K. C., Blundell K. M., 2015, A&A, 574, A30
  * Roberts et al. (2018) Roberts D. H., Saripalli L., Wang K. X., Sathyanarayana Rao M., Subrahmanyan R., KleinStern C. C., Morii-Sciolla C. Y., Simpson L., 2018, ApJ, 852, 47
  * Romanova & Lovelace (1992) Romanova M. M., Lovelace R. V. E., 1992, A&A, 262, 26
  * Ruszkowski et al. (2007) Ruszkowski M., Enßlin T. A., Brüggen M., Heinz S., Pfrommer C., 2007, MNRAS, 378, 662
  * Rybicki & Lightman (1979) Rybicki G. B., Lightman A. P., 1979. Wiley, New York
  * Ryu et al. (2000) Ryu D., Jones T. W., Frank A., 2000, ApJ, 545, 475
  * Savini et al. (2018) Savini F., et al., 2018, MNRAS, 474, 5023
  * Scheuer (1974) Scheuer P. A. G., 1974, MNRAS, 166, 513
  * Scheuer (1982) Scheuer P. A. G., 1982, Extragalactic radio sources; IAU Symposium, 97, 163
  * Shulevski et al. (2017) Shulevski A., et al., 2017, A&A, 600, A65
  * Sikora et al. (2005) Sikora M., Begelman M. C., Madejski G. M., Lasota J.-P., 2005, ApJ, 625, 72
  * Sironi & Spitkovsky (2014) Sironi L., Spitkovsky A., 2014, ApJ, 783, L21
  * Sironi et al. (2015) Sironi L., Petropoulou M., Giannios D., 2015, MNRAS, 450, 183
  * Smith & Donohoe (2019) Smith M. D., Donohoe J., 2019, MNRAS, 490, 1363
  * Smith & Donohoe (2021) Smith M. D., Donohoe J., 2021, MNRAS, 502, 423
  * Steenbrugge et al. (2010) Steenbrugge K. C., Heywood I., Blundell K. M., 2010, MNRAS, 401, 67
  * Tchekhovskoy (2015) Tchekhovskoy A., 2015, in Contopoulos I., Gabuzda D., Kylafis N., eds, Astrophysics and Space Science Library Vol. 414, The Formation and Disruption of Black Hole Jets. p. 45, doi:10.1007/978-3-319-10356-3_3
  * Tout & Pringle (1996) Tout C. A., Pringle J. E., 1996, MNRAS, 281, 219
  * Tremblay et al. (2010) Tremblay G. R., O’Dea C. P., Baum S. A., Koekemoer A. M., Sparks W. B., de Bruyn G., Schoenmakers A. P., 2010, ApJ, 715, 172
  * Turk et al. (2011) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, The Astrophysical Journal Supplement Series, 192, 9
  * Vaidya et al. (2016) Vaidya B., Mignone A., Bodo G., Massaglia S., 2016, in Journal of Physics Conference Series. p. 012023, doi:10.1088/1742-6596/719/1/012023
  * Vaidya et al. (2018) Vaidya B., Mignone A., Bodo G., Rossi P., Massaglia S., 2018, ApJ, 865, 144
  * Vernaleo & Reynolds (2006) Vernaleo J. C., Reynolds C. S., 2006, ApJ, 645, 83
  * Weinberger et al. (2017) Weinberger R., Ehlert K., Pfrommer C., Pakmor R., Springel V., 2017, MNRAS, 470, 4530
  * Werner et al. (2010) Werner N., et al., 2010, MNRAS, 407, 2063
  * Yates-Jones et al. (2021) Yates-Jones P. M., Shabala S. S., Krause M. G. H., 2021, MNRAS, 508, 5239
  * Young et al. (2005) Young A. J., Wilson A. S., Tingay S. J., Heinz S., 2005, ApJ, 622, 830
  * Zhuravleva et al. (2015) Zhuravleva I., et al., 2015, MNRAS, 450, 4184