# Blandford-Znajek jets in galaxy formation simulations: method and
implementation

Rosie Y. Talbot1, Martin A. Bourne1,2 and Debora Sijacki1,2
1Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge,
CB3 0HA, UK
2Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge, CB3 0HA, UK E-mail<EMAIL_ADDRESS>(RYT)

(Submitted to MNRAS)

###### Abstract

Jets launched by active galactic nuclei (AGN) are believed to play a
significant role in shaping the properties of galaxies and provide an
energetically viable mechanism through which galaxies can become quenched.
Here we present a novel AGN feedback model, which we have incorporated into
the arepo code, that evolves the black hole mass and spin as the accretion
flow proceeds through a thin $\alpha$-disc which we self-consistently couple
to a Blandford-Znajek jet. We apply our model to the central region of a
typical radio-loud Seyfert galaxy embedded in a hot circumgalactic medium
(CGM). We find that jets launched into high pressure environments thermalise
efficiently due to the formation of recollimation shocks and the vigorous
instabilities that these shocks excite increase the efficiency of the mixing
of CGM and jet material. The beams of more overpressured jets, however, are
not as readily disrupted by instabilities so the majority of the momentum flux
at the jet base is retained out to the head, where the jet terminates in a
reverse shock. All jets entrain a significant amount of cold circumnuclear
disc material which, while energetically insignificant, dominates the lobe
mass together with the hot, entrained CGM material. The jet power evolves
significantly due to effective self-regulation by the black hole, fed by
secularly-driven, intermittent mass flows. The direction of jets launched
directly into the circumnuclear disc changes considerably due to effective
Bardeen-Petterson torquing. Interestingly, these jets obliterate the innermost
regions of the disc and drive large-scale, multi-phase, turbulent, bipolar
outflows.

###### keywords:

black hole physics – methods: numerical – galaxies: jets – galaxies: active

††pubyear: 2020††pagerange: Blandford-Znajek jets in galaxy formation
simulations: method and implementation–B

##  1 Introduction

In the centres of many galaxy clusters the hot, X-ray emitting gas is observed
to have rapid ($\lesssim 1\;{\rm Gyr}$) cooling times. Without invoking some
heating mechanism, significant quantities of cold gas should flow towards the
centre of these clusters (Peterson & Fabian, 2006; Fabian, 1994). Many such
clusters, however, do not show evidence of high levels of star formation
(Cooke et al., 2016; Donahue et al., 2010) nor do they have large reservoirs
of molecular gas (Fogarty et al., 2019; Russell et al., 2019, 2017), both of
which should be present if such a cooling flow were to exist. It naturally
follows that some sort of heating process must be at play, offsetting
radiative cooling and ensuring that the intracluster medium (ICM) remains in
approximate thermal equilibrium.

Jet feedback from the AGN associated with the central brightest cluster galaxy
is often invoked as such a mechanism (McNamara & Nulsen, 2012, 2007) and in
many cases there is direct observational evidence of the interaction of jet
lobes with the surrounding medium with estimated cavity energies comparable to
that required to offset catastrophic cooling. Several alternative physical
mechanisms by which energy can be imparted to the ICM have been proposed such
as thermal conduction (Ruszkowski & Oh, 2010; Voigt & Fabian, 2004; Zakamska &
Narayan, 2003; Narayan & Medvedev, 2001). In reality the heating of the ICM is
likely due to the collective action of a multitude of processes although their
relative importance remains to be determined.

AGN jet feedback, however, is not confined just to these most massive of
systems. Indeed, there is growing evidence, both observational and
theoretical, that AGN feedback is acting in a variety of galaxy contexts,
ranging from spirals (e.g. Mao et al., 2015; Croston et al., 2008; Ledlow et
al., 2001) to dwarfs (e.g. Koudmani et al., 2019; Mezcua et al., 2019; Mezcua
et al., 2018; Chilingarian et al., 2018; Penny et al., 2018; Silk, 2017).
Galaxy formation simulations, however, typically disfavour strong AGN feedback
in low mass galaxies and, specifically, the role of AGN jets in dwarfs and
spirals is largely unexplored.

There is clear observational evidence for radio jets in Seyfert galaxies
(Williams et al., 2017; Mingo et al., 2012, 2011; Gallimore et al., 2006; Hota
& Saikia, 2006; Morganti et al., 1999). Interestingly, some studies find
little evidence for correlation between the jet axis and the major axis of the
host galaxy in these systems (Nagar & Wilson, 1999; Clarke et al., 1998;
Schmitt et al., 1997) and some show signs of ‘bent’ radio jets which could be
due to jet precession or interaction with their surroundings. The lack of
correlation could come about due to the large-scale gas angular momentum not
being aligned with that of the innermost regions of the accretion disc (Kinney
et al., 2000; Pringle, 1997, 1996). Upcoming new observational data from Lynx,
Athena and SKA will shed new light on these issues by allowing us to
investigate these energetically moderate sources in much greater detail both
locally and at high redshifts.

There are clearly numerous open questions surrounding the nature of radio-mode
feedback: ‘Can the energy from the jets be efficiently and (largely)
isotropically communicated to the surrounding medium and if so, what are the
dominant physical mechanisms that are responsible?’, ‘How are these jets
launched and what is their duty cycle?’, ‘What processes drive the transition
to quasar-mode feedback?’ to name but a few.

A large body of work studying X-ray black hole binaries has provided crucial
insights as to the transition between low and high accretion state (Begelman &
Armitage, 2014; Fender et al., 2004). Theoretical models of AGN accretion and
feedback that aim to motivate the transition from radio to quasar mode
feedback have been developed in analogy with this (see e.g. Yuan & Narayan,
2014; Merloni & Heinz, 2008; Churazov et al., 2005). Direct observational
evidence, however, is much less conclusive due to the significantly longer
timescales associated with supermassive black hole accretion relative to those
of X-ray binary systems.

Simulations of AGN jets prove a vital tool to shed further light on these
questions. The processes expected to be at play, however, span a vast range of
scales and, hence, many different ‘flavours’ of simulations have developed,
each of which is specialised to probe the relevant physics in different
regimes.

On the smallest scales, general relativistic magnetohydrodynamic (GRMHD)
simulations have proved an effective tool to investigate the formation of
these jets and how they interact with the accretion flow that feeds the black
hole (Chatterjee et al., 2019; Liska et al., 2018; Sądowski & Narayan, 2015;
Tchekhovskoy et al., 2011; Penna et al., 2010; McKinney, 2006). In these
simulations, crucially, the jet launching scale is resolvable and some
simulations find that evacuated funnels of outflowing magnetically dominated
plasma develop self-consistently along the black hole spin axis. Due to
computational constraints, these simulations mainly focus on thick discs which
are expected to form when the accretion rate is significantly sub or super-
Eddington. It is only very recently, thanks to significant computational
advances, that the thin-disc paradigm has been started to be explored. For
example, Liska et al. (2020, 2019) have shown that the inner regions of
initially misaligned, thin discs are torqued into alignment with the
equatorial plane of the black hole, leading to the formation of warped discs
(in accordance with the predictions of Bardeen & Petterson (1975)). The fact
that these thin discs are able to launch precessing, relativistic, magnetised
jets raises the question of whether thick discs are a necessary requirement
for the launching of radio jets.

Moving up to larger scales, there are numerous studies that examine jet
launching in hydrostatic atmospheres and idealised cluster setups, with the
aim of better understanding the jet-ICM interaction (see e.g. Bambic &
Reynolds, 2019; Ehlert et al., 2018; Bourne & Sijacki, 2017; Weinberger et
al., 2017b; Yang & Reynolds, 2016b; Reynolds et al., 2015; Krause et al.,
2012; Omma et al., 2004). Due to the large dynamic range needed to capture the
disparate scales associated with this problem, it is not possible to entirely
self-consistently launch jets; some sub-grid prescription must be used
instead. Many of these studies fix the jet power and direction, which does not
allow for the formation of self-regulated feedback loop, making it hard to
reconcile these results with the findings of GRMHD simulations.

AGN feedback in galaxy formation simulations has been widely considered (Davé
et al., 2019; Henden et al., 2018; Weinberger et al., 2018; McCarthy et al.,
2017; Sijacki et al., 2015; Schaye et al., 2015; Crain et al., 2015; Dubois et
al., 2012; Sijacki et al., 2007). The inclusion of supermassive black holes in
these simulations has proved necessary in order to regulate the properties of
massive galaxies, thus producing more realistic galaxy populations that are
broadly consistent with observations. In these simulations, however, the
typical resolution scale is usually on the order of $100\;{\rm pc}$ at best
and so it is often difficult to resolve the lobe inflation process, let alone
the jet launching process. This usually leads to the jet energy being injected
in the form of a thermal or kinetic dump, e.g. in off-centre bubbles or large-
scale (bipolar) winds, which eliminates the requirement that the lobe
inflation be followed, albeit at the cost of reduced physical accuracy. Zoom-
in simulations of cosmological clusters somewhat address the resolution
problems associated with full cosmological box simulations. Indeed, individual
jet studies that follow the lobe inflation process in a live cosmological
environment have been carried out (Bourne & Sijacki, 2020; Bourne et al.,
2019; Dubois et al., 2010; Morsony et al., 2010; Heinz et al., 2006).

Despite these problems, the modelling of jet feedback in galaxy and
cosmological scale simulations is becoming more sophisticated, incorporating
more accurate modelling of (computationally) complex physical processes such
as magnetic fields and cosmic rays (Yang et al., 2019; Ehlert et al., 2018;
Weinberger et al., 2017b). Jet studies are also now finding that they are able
to reproduce features seen in radio and X-ray observations with unprecedented
accuracy (Bourne & Sijacki, 2020; Bourne et al., 2019). With these advances in
mind it is timely to ensure that self-regulation of the jet feedback is fully
investigated in these contexts. Whilst including self-regulated jet feedback
in galaxy-scale simulations is by no means a new idea (Yang & Reynolds, 2016b;
Li & Bryan, 2014b; Gaspari et al., 2011; Cattaneo & Teyssier, 2007; Vernaleo &
Reynolds, 2006), self-consistent models that link the evolution of gas angular
momentum with the jet power are still in their infancy.

Due to the vast range of scales that are involved in tracking accretion flows
and AGN jet propagation in galaxy-scale simulations, many of the launching
processes (that can be captured to some extent in GRMHD simulations) fall
below scales that are resolvable. Our work here focuses on galaxy-scale
simulations but aims to (somewhat) bridge the gap between galaxy-scale and
GRMHD simulations. Since we cannot resolve the jet launching scales we have
developed a self-consistent sub-grid model for AGN feedback in the form of a
Blandford-Znajek jet, motivated by the results of GRMHD simulations and
general relativistic analytic predictions. Our model assumes that the black
hole is surrounded by a sub-grid, thin $\alpha$-disc which modulates the
accretion flow onto the black hole, allowing us to accurately follow both the
evolution of the mass and angular momentum of the black hole. In such a model
it is of vital importance to self-consistently evolve the system as a whole
and hence the evolution of the black hole as a result of launching the jet is
included in our model as well.

Our model has been designed such that it can be employed in galaxy scale
simulations with the ultimate aim of assessing whether self-consistent
Blandford-Znajek jets can reproduce large-scale galaxy and cluster
observables. This work, however, specifically focuses on resolving the parsec-
scale interactions of jets with the cold dense circumnuclear gas found close
to the black hole as well as the interactions of the jet with the surrounding
hot CGM.

The structure of this paper is as follows. In Section 2 we motivate and
develop the analytic equations that govern the evolution of our sub-grid
model. Then in Section 3 we describe how this is implemented into arepo and
explain the refinement schemes we employ. In Section 4 we introduce the
simulation setup that we use to test our sub-grid model. In Section 5 we
discuss simulations in which the jet power and direction are fixed which we
then use to better understand the results of analogous simulations in which we
utilise our full Blandford-Znajek jet model, which are presented in Section 6.
In Section 7 we discuss our work in the context of previous galaxy-scale
simulations and highlight any physics missing from our simulations and
finally, we end with our conclusions in Section 8.

##  2 Theoretical background

In our sub-grid prescription for black hole accretion and feedback, we model
the black hole as a sink particle surrounded by a sub-grid accretion disc
which follows the structure of a steady-state, geometrically-thin
$\alpha$-disc, as presented in Shakura & Sunyaev (1973). The black hole has a
sub-resolution mass of $M_{\rm BH}$ and angular momentum $\mathbf{J}_{\rm
BH}$. Similarly, the disc is characterised solely by its mass, $M_{\rm d}$,
and total angular momentum, $\mathbf{J}_{\rm d}$. The sub-grid disc-black hole
system interacts gravitationally with the wider simulation via its dynamical
mass, $M_{\rm dyn}=M_{\rm BH}+M_{\rm d}$.

In order to follow the evolution of the system and the launching of the jet,
we need to determine equations that specify the evolution of the mass and
angular momentum of the black hole and disc as well as equations that
determine the mass, energy and momentum fluxes into the jet. We break down
this problem by considering the physical processes that will lead to evolution
of the system. Specifically, in our model, the mass-energy of the black hole
evolves due to accretion from the sub-grid disc and the energy loss associated
with launching the Blandford-Znajek jet. The mass of the disc evolves due to
inflow from the wider simulation and mass that leaves the disc at its inner
edge. The angular momentum of the black hole evolves due to accretion from the
sub-grid disc, the Bardeen-Petterson effect and the launching of the
Blandford-Znajek jet. Finally, the angular momentum of the disc evolves due to
inflow from the wider simulation, the mass that leaves its inner edge and the
Bardeen-Petterson effect.

Our implementation of accretion via a thin, potentially warped $\alpha$-disc
largely follows the procedure presented in Fiacconi et al. (2018). Our
numerical implementation, however, is somewhat different (see Section 3) and
the launching of the Blandford-Znajek jet necessarily changes the equations
governing the evolution of the system. We now proceed to fully develop these
evolutionary equations.

###  2.1 Accretion and inflow

We first determine how the black hole and disc will evolve due to inflow of
material from the wider simulation onto the sub-grid disc and the flow of mass
onto the black hole from the inner edge of the disc.

In the Shakura & Sunyaev (1973) model, the disc extends down the the innermost
stable circular orbit (ISCO) of the black hole and accretion is driven by an
effective viscosity, $\nu$, parameterised by $\alpha$, according to

$\nu\approx\alpha c_{\rm s}\,H\,,$ (1)

where $c_{\rm s}$ is the sound speed and $H$ is the disc scale-height. Whilst
the value of this $\alpha$-parameter is poorly constrained by theory, it
enters into the structural equations raised only to low powers. We therefore
choose to fix it at $0.1$, consistent with some observations (King et al.,
2007; Schreiber et al., 2003; Cannizzo, 2001a, b).

The mass of the black hole evolves due to the accretion of
material111Specifically the binding energy associated with the material. from
the inner edge of the $\alpha$-disc. If the rest-mass flux onto the black hole
is $\dot{M}_{\rm{BH},0}$ then the corresponding increase in black hole mass
will be

$\dot{M}_{\rm BH}^{\rm(acc)}=(1-\epsilon_{\rm r})\,\dot{M}_{\rm BH,0}\,,$ (2)

where $\epsilon_{\rm r}$ is the spin dependent radiative efficiency which
corrects for the reduced binding energy of material passing through the ISCO
(Novikov & Thorne, 1973) (see Appendix A for the full expression for
$\epsilon_{\rm r}$).

The mass of the accretion disc will evolve due to inflow of material from the
surroundings as well as the loss of material at its inner edge. Assuming that
the inflowing material circularises and settles into the disc, $M_{\rm d}$
will change according to

$\dot{M}_{\rm d}^{\rm(acc)}=\dot{M}_{\rm in}-\dot{M}\,,$ (3)

where $\dot{M}_{\rm in}$ is the mass inflow rate from the surroundings onto
the black hole-disc system and $\dot{M}$ is the steady state rest mass flux
through the disc.

In the above equations we have taken care to differentiate between the rest
mass flux onto the black hole, $\dot{M}_{\rm BH,0}$ and the mass flux through
the disc, $\dot{M}$. This allows for the possibility that not all material
that leaves the disc at its inner edge is captured by the black hole. We do
this as later we will assume that the jet is mass loaded with material from
the accretion disc (see Section 2.4).

The angular momentum direction of the accreted material is fixed by the
assumption that the Bardeen-Petterson effect causes the inner disc to align
with the equatorial plane of the black hole (as discussed in Section 2.2
below). This means that the angular momentum evolution of the black hole due
to these mass flows is

$\dot{\mathbf{J}}_{\rm BH}^{\rm(acc)}=L_{\rm ISCO}\;\dot{M}_{\rm
BH,0}\;\mathbf{j}_{\rm BH}\,,$ (4)

where $L_{\rm ISCO}$ is the magnitude of the specific angular momentum of the
gas that is accreted by the black hole from the disc (see Appendix A for the
full expression of $L_{\rm ISCO}$ in terms of the black hole spin) and
$\mathbf{j}_{\rm BH}$ is the unit vector in the direction of the black hole
angular momentum.

Similarly, the disc angular momentum evolution is given by

$\dot{\mathbf{J}}_{\rm d}^{\rm(acc)}=\dot{M}_{\rm in}\mathbf{L}_{\rm
in}-L_{\rm ISCO}\dot{M}\,\mathbf{j}_{\rm BH}\;,$ (5)

where $\mathbf{L}_{\rm in}$ is the specific angular momentum of the inflowing
gas.

In accordance with the assumptions of the Shakura & Sunyaev (1973) thin disc
equations we cap the mass flux through the disc, $\dot{M}$, at the Eddington
rate, $\dot{M}_{\rm Edd}$. Parameterising $\dot{M}$ in terms a fraction
$f_{\rm Edd}$ of $\dot{M}_{\rm Edd}$ we can write

$\displaystyle\dot{M}$ $\displaystyle={\rm min}(f_{\rm Edd},1)\;\dot{M}_{\rm
Edd}\,,$ $\displaystyle={\rm min}(f_{\rm Edd},1)\;\frac{M_{\rm BH}}{\tau_{\rm
Salp}}\;,$ (6)

where $\tau_{\rm Salp}$ is the Salpeter time

$\tau_{\rm Salp}=\frac{\epsilon_{\rm r}\,\kappa_{\rm ES}c}{4\pi G}\approx
4.5\times 10^{7}\,\bigg{(}\frac{\epsilon_{\rm r}}{0.1}\bigg{)}\;{\rm yrs}\;,$
(7)

and $\kappa_{\rm ES}\approx 0.4\;{\rm cm^{2}g^{-1}}$ is the electron-
scattering opacity. $f_{\rm Edd}$ is determined using the solutions to the
Shakura & Sunyaev (1973) thin disc equations under the assumption that the
disc is gas pressure dominated, in steady state and that electron-scattering
dominates the opacity which altogether give

$f_{\rm Edd}\approx 0.76\bigg{(}\frac{\epsilon_{\rm
r}}{0.1}\bigg{)}\bigg{(}\frac{M_{\rm d}}{10^{4}\,{\rm
M_{\odot}}}\bigg{)}^{5}\bigg{(}\frac{M_{\rm BH}}{10^{6}\,{\rm
M_{\odot}}}\bigg{)}^{-47/7}\bigg{(}\frac{a\;J_{\rm d}/J_{\rm
BH}}{3}\bigg{)}^{-25/7}\,,$ (8)

where $a$ is the dimensionless spin parameter

$a=\frac{c\,J_{\rm BH}}{G\,M_{\rm BH}^{2}}\,.$ (9)

Theory constrains $a$ to lie in the range $0\leq a\leq 1$, however, if the
black hole is surrounded by a thin disc then preferential capture of negative
angular momentum photons from the surface of the disc will cap the spin such
that $a\leq 0.998$ (Thorne, 1974).

At low accretion rates we expect radiative cooling to become inefficient and
the disc to thicken (Narayan & Yi, 1994, 1995; Abramowicz et al., 1995) and in
this regime the Shakura & Sunyaev (1973) model breaks down. We therefore
ensure that the accretion rate does not become too low by imposing $f_{\rm
Edd}>f_{\rm Edd}^{\rm(min)}=10^{-4}$, as in Fiacconi et al. (2018).

Equation (3) admits the possibility that the accretion disc could grow without
bound if $\dot{M}_{\rm in}$ is consistently higher than $\dot{M}$. Were this
to be the case, the disc could enter a state in which self-gravity becomes
important in its outer regions. This self-gravitating regime is highly non-
linear and could lead to processes such as the formation of spiral arms and
bars, fragmentation of the disc and, ultimately, star formation. This regime
is poorly understood (Goodman, 2003) and it is not known how these processes
may couple to alter the mass and angular momentum transport through the disc.
We therefore choose to prevent the disc from entering such a regime by capping
the mass inflow rate such that the disc mass remains below the mass at which
the material at the outer edge of the disc would be unstable according to the
Toomre criterion

$M_{\rm d}^{\rm(sg)}\approx 2.2\times 10^{4}\;\bigg{(}\frac{f_{\rm
Edd}}{\epsilon_{\rm r}/0.1}\bigg{)}^{4/45}\bigg{(}\frac{M_{\rm
BH}}{10^{6}\,{\rm M_{\odot}}}\bigg{)}^{34/45}\;{\rm M_{\odot}}\,.$ (10)

For the derivation of the above expression for $M_{\rm d}^{\rm(sg)}$, see
appendix A of Fiacconi et al. (2018).

###  2.2 The Bardeen-Petterson effect

A spinning black hole leads to gas on non-equatorial orbits precessing about
the angular momentum vector of the black hole at a rate that depends on the
distance from the black hole, i.e. the Lense-Thirring effect (Lense &
Thirring, 1918). This precession is countered by the vertical viscosity in the
disc which, in the thin disc regime, can then lead to the disc becoming
warped. The inner regions of an initially misaligned disc will (counter-)align
with the spin of the black hole while the outer disc remains misaligned, the
so-called Bardeen-Petterson effect (Bardeen & Petterson, 1975). The radius at
which the transition occurs, the warp radius $R_{\rm warp}$, is determined by
a balance of the vertical warp propagation timescale, $\tau_{\nu_{2}}$, and
the Lense-Thirring precession timescale, $\tau_{\rm LT}$,

$\frac{R_{\rm warp}^{2}}{\nu_{2}}\equiv\tau_{\nu_{2}}(R_{\rm warp})=\tau_{\rm
LT}(R_{\rm warp})\equiv\frac{c^{2}R_{\rm warp}^{3}}{2\,G\,J_{\rm BH}}\,,$ (11)

where $\nu_{2}$ is the vertical viscosity (see e.g. Dotti et al., 2013; Martin
et al., 2007; Lodato & Pringle, 2006). In the steady, thin disc regime the
warp will propagate outwards diffusively and, provided the warp radius lies
within the disc, the disc can attain a steady, warped state since the warping
timescale is shorter than the local viscous timescale. Gas in the misaligned
outer disc will experience a torque from the black hole that is strongest in
the region around the warp radius. In an isolated black hole-disc system, the
total angular momentum, $\mathbf{J}_{\rm tot}=\mathbf{J}_{\rm
BH}+\mathbf{J}_{\rm d}$ is conserved, meaning that an equal and opposite
torque must be felt by the black hole itself in order that $\mathbf{J}_{\rm
tot}$ is conserved. Altogether, this leads to precession and
(counter-)alignment of the black hole and disc angular momentum vectors.

Since the torque exerted on the black hole due to the warping of the disc does
not have a component in the direction of the black hole angular momentum
(Pringle, 1992) it can be expressed in the form

$\dot{\mathbf{J}}_{\rm BH}^{\rm(BP)}=-J_{\rm
BH}\bigg{\\{}K_{1}(\mathbf{j}_{\rm BH}\times\mathbf{j}_{\rm
d})+K_{2}\big{[}\mathbf{j}_{\rm BH}\times(\mathbf{j}_{\rm
BH}\times\mathbf{j}_{\rm d})\big{]}\bigg{\\}}\;,$ (12)

where $\mathbf{j}_{\rm BH}$ and $\mathbf{j}_{\rm d}$ are the unit vectors in
the direction of the black hole and disc angular momenta, respectively. The
first term in equation (12) corresponds to a precession about $\mathbf{J}_{\rm
tot}$ and the second term leads to alignment with $\mathbf{J}_{\rm tot}$ (King
et al., 2005). In general $K_{1}$ and $K_{2}$ are arbitrary functions,
however, under the assumptions that (i) the misalignment between the angular
momenta is small and (ii) the total angular momentum is dominated by that of
the disc, we can use the results of Martin et al. (2007) to identify $K_{1}$
and $K_{2}$ with the precession and alignment timescales. Whilst these
assumptions are inherently restrictive, we find that this is the best option
available to determine the Bardeen-Petterson torque, short of solving for the
full structure of the accretion disc. We discuss the implications of this
assumption in Section 6.

Ultimately, this allows us to express $K_{1}$ and $K_{2}$ as

$K_{1}=\frac{\sin(\pi/7)}{\tau_{\rm GM}}\qquad\text{and}\qquad
K_{2}=\frac{\cos(\pi/7)}{\tau_{\rm GM}}\,,$ (13)

where $\tau_{\rm GM}$ is the gravitomagnetic timescale:

$\tau_{\rm GM}\approx 0.17\bigg{(}\frac{M_{\rm BH}}{10^{6}\,{\rm
M_{\odot}}}\bigg{)}^{-2/35}\bigg{(}\frac{f_{\rm Edd}}{\epsilon_{\rm
r}/0.1}\bigg{)}^{-32/35}a^{5/7}\;\rm{Myrs}\;,$ (14)

(Dotti et al., 2013; Perego et al., 2009; Martin et al., 2007).

Equation (12) and the identities in equation (13) together determine the
evolution of the black hole angular momentum due to the Bardeen Petterson
effect. The evolution of the disc is then simply

$\dot{\mathbf{J}}_{\rm d}^{\rm(BP)}=-\dot{\mathbf{J}}_{\rm BH}^{\rm(BP)}\;.$
(15)

If the warp radius is larger than the disc extent, however, the disc will not
reach a steady warped state and these equations for the angular momentum
evolution are no longer valid. This will be the case if the mass of the black
hole exceeds a characteristic warping mass222For the derivation of this
warping mass see appendix A of Fiacconi et al. (2018).

$M_{\rm BH}^{\rm(warp)}\approx 10^{7}\,\bigg{(}\frac{M_{\rm d}}{10^{4}\,{\rm
M_{\odot}}}\bigg{)}^{35/82}\,\bigg{(}\frac{f_{\rm Edd}}{\epsilon_{\rm
r}/0.1}\bigg{)}^{-17/82}\;a^{-25/82}\;{\rm M_{\odot}}\;.$ (16)

In this scenario the disc angular momentum will (counter-)align with that of
the black hole on a timescale set by the diffusive warp propagation. In this
regime we choose to employ a simple prescription and assume that the disc
angular momentum instantaneously (counter-)aligns with that of the black hole.
Alignment will occur if $\mathbf{j}_{\rm BH}\cdot\mathbf{j}_{\rm d}>-J_{\rm
BH}/2J_{\rm d}$ (King et al., 2005), otherwise they will counter-align.

###  2.3 The Blandford-Znajek mechanism

Blandford & Znajek (1977) proposed a mechanism through which magnetic fields
could extract rotational energy from the black hole. A large poloidal magnetic
field that threads a spinning black hole becomes twisted due to frame-dragging
effects and a toroidal field develops. The resulting magnetic pressure in this
configuration then provides the launching mechanism for the jet with the work
done by the black hole on the field lines leading to the extraction of its
rotational energy.

Blandford & Znajek (1977) determined the perturbations to the fields under the
assumption that the black hole is surrounded by a razor-thin disc and its
magnetosphere is force-free. If the magnetic flux threading one hemisphere of
the black hole horizon is

$\Phi_{\rm BH}=\frac{1}{2}\int\left|B^{r}\right|\mathrm{d}A_{\theta\phi}\,,$
(17)

where $B^{r}$ is the radial component of the magnetic field and
$\mathrm{d}A_{\theta\phi}$ is the area element on the horizon then, in the
limit of $a\ll 1$ and with fixed $\Phi_{\rm BH}$, the rotational energy of the
black hole can be extracted at a rate

$\dot{E}_{\rm BZ}=\frac{\kappa}{4\pi}\Phi_{\rm BH}^{2}\frac{a^{2}\,c}{16r_{\rm
g}^{2}}\,,$ (18)

where $r_{\rm g}=G\,M_{\rm BH}/c^{2}$ is the black hole gravitational radius
and $\kappa$ is a parameter that depends only weakly on the magnetic field
geometry. In this work we fix $\kappa=1/(6\pi)$ which corresponds to a split-
monopole geometry.

Launching a Blandford-Znajek jet leads to a decrease in the mass-energy of the
black hole and causes it to spin-down. The energy flux due to the Blandford-
Znajek effect has been expanded analytically to higher orders under various
assumptions (Tchekhovskoy et al., 2010; Tanabe & Nagataki, 2008) and there
have been many studies investigating the process via GRMHD simulations (Penna
et al., 2013; McKinney et al., 2012; Hawley & Krolik, 2006; McKinney & Gammie,
2004; Komissarov, 2001). In Appendix B we provide an explicit derivation of
the energy and angular momentum fluxes to second and first order in the
perturbed fields, respectively. When implementing our sub-grid model, however,
we also include the higher order terms in the energy flux determined
numerically by Tchekhovskoy et al. (2010) as this improves the accuracy of the
determination of the energy flux at higher spin values. We can therefore write
the energy flux as

$\dot{E}_{\rm BZ}=\frac{\kappa}{4\pi}\,\phi_{\rm BH}^{2}\,f(a)\,\dot{M}_{\rm
BH,0}\,c^{2}\,,$ (19)

where the dimensionless magnetic flux is

$\phi_{\rm BH}=\frac{\Phi_{\rm BH}}{\sqrt{\dot{M}_{\rm BH,0}\,r_{\rm
g}^{2}c}}\,,$ (20)

and $f(a)$ is a dimensionless function of the black hole spin, given by

$\displaystyle
f(a)=\Bigg{(}\frac{a}{2\,\big{(}1+\sqrt{1-a^{2}}\big{)}}\Bigg{)}^{2}+1.38\,$
$\displaystyle\Bigg{(}\frac{a}{2\,\big{(}1+\sqrt{1-a^{2}}\big{)}}\Bigg{)}^{4}$
$\displaystyle-9.2\,$
$\displaystyle\Bigg{(}\frac{a}{2\,\big{(}1+\sqrt{1-a^{2}}\big{)}}\Bigg{)}^{6}\,.$
(21)

This outward flux of energy leads to the evolution of black hole mass
according to

$\dot{M}_{\rm BH}^{\rm(BZ)}=-\epsilon_{\rm BZ}\dot{M}_{\rm BH,0}c^{2}\,,$ (22)

where, in analogy with the radiative efficiency, we have introduced a
‘Blandford-Znajek efficiency’

$\epsilon_{\rm BZ}\equiv\frac{\kappa}{4\pi}\phi_{\rm BH}^{2}f(a)\,.$ (23)

The corresponding outward angular momentum flux leads to evolution of the
black hole angular momentum

$\dot{\mathbf{J}}_{\rm BH}^{\rm(BZ)}=-\dot{M}_{\rm BH,0}\,\mathbf{L}_{\rm
BZ}\,,$ (24)

which we have defined in terms of an effective specific angular momentum

$\displaystyle\mathbf{L}_{\rm BZ}$ $\displaystyle=L_{\rm BZ}\,\mathbf{j}_{\rm
BH}\;,$ $\displaystyle=\frac{\kappa}{2\pi}\,\phi_{\rm
BH}^{2}\,\Bigg{(}\frac{a}{2\,\big{(}1+\sqrt{1-a^{2}}\big{)}}\Bigg{)}\,\frac{G\,M_{\rm
BH}}{c}\,\mathbf{j}_{\rm BH}\,.$ (25)

The dimensionless magnetic flux that threads the black hole horizon,
$\phi_{\rm BH}$, cannot be directly measured from our simulations since we do
not track magnetic fields (and even if we did, the processes that determine it
would be well below our resolution scale). We therefore choose to use the spin
dependent values for $\phi_{\rm BH}$ obtained from the GRMHD simulations in
Tchekhovskoy et al. (2012). It should be noted here that these simulations
correspond to an accretion disc in a magnetically arrested state (i.e. the
magnetic flux on the black hole has reached saturation). We will discuss this
choice further in Section 7.

###  2.4 Jet launching

When launching the Blandford-Znajek jet we explicitly mass load it by assuming
that some of the mass that flows through the ISCO of the sub-grid accretion
disc is drawn up into the jet. We control the fraction of this mass flux that
goes into the jet using the mass-loading factor, $\eta_{\rm J}$ (Bourne &
Sijacki, 2017; Ostriker et al., 2010), a free parameter, which we define to be

$\dot{M}_{\rm J}=\frac{\eta_{\rm J}}{1+\eta_{\rm J}}\dot{M}\,.$ (26)

Throughout this work we choose $\eta_{\rm J}=1$ (Bourne & Sijacki, 2017).

In terms of the mass-loading factor, we can express the rest mass flux across
the black hole horizon as

$\dot{M}_{\rm BH,0}=\frac{1}{1+\eta_{\rm J}}\dot{M}\,.$ (27)

The energy and momentum fluxes into the jet are

$\displaystyle\dot{E}_{\rm J}$ $\displaystyle=\frac{1}{2}\dot{M}_{\rm J}v_{\rm
J}^{2}\,,$ (28) $\displaystyle\dot{P}_{\rm J}$ $\displaystyle=\dot{M}_{\rm
J}v_{\rm J}\,,$ (29)

where $v_{\rm J}$ is the sub-resolution velocity of the jet.

We assume that the jet is powered by the entire outward flux of energy from
the black hole, as predicted by the Blandford-Znajek process, and thus

$\displaystyle\dot{E}_{\rm J}$ $\displaystyle=\epsilon_{\rm BZ}\dot{M}_{\rm
BH,0}c^{2}\,,$ (30) $\displaystyle\dot{P}_{\rm J}$
$\displaystyle=\sqrt{2\dot{E}_{\rm J}\dot{M}_{\rm J}}=\frac{\sqrt{2\eta_{\rm
J}\,\epsilon_{\rm BZ}}}{1+\eta_{\rm J}}\dot{M}c\,.$ (31)

It should be noted here that in our numerical scheme we inject the jet energy
into a finite mass. Since this leads to numerical mass-loading, the energy and
momentum fluxes into the jet cannot satisfy equations (30) and (31)
simultaneously. We choose to conserve kinetic energy and so the momentum flux
in equation (31) will not correspond to the momentum flux into the jet and
similarly the sub-resolution velocity, $v_{\rm J}$, will not correspond to the
initial velocity of jet material in the simulations. The specifics of the jet
injection are detailed in Section 3.6 below.

###  2.5 Summary of equations

We end this section by collating the relevant equations, developed above, into
a set of overarching equations that govern the evolution of the black hole and
accretion disc and the properties of the jet in our model.

The black hole mass evolves according to

$\displaystyle\dot{M}_{\rm BH}$ $\displaystyle=\dot{M}_{\rm
BH}^{\rm(acc)}+\dot{M}_{\rm BH}^{\rm(BZ)}\,,$ (32)
$\displaystyle=\frac{(1-\epsilon_{\rm r}-\epsilon_{\rm BZ})}{1+\eta_{\rm
J}}\dot{M}\,,$ (33)

and the disc mass

$\displaystyle\dot{M}_{\rm d}=\dot{M}_{\rm in}-\dot{M}\,.$ (34)

The black hole angular momentum is given by

$\displaystyle\dot{\mathbf{J}}_{\rm BH}=$
$\displaystyle\,\dot{\mathbf{J}}_{\rm BH}^{\rm(acc)}+\dot{\mathbf{J}}_{\rm
BH}^{\rm(BZ)}+\dot{\mathbf{J}}_{\rm BH}^{\rm(BP)}\,,$ (35) $\displaystyle=$
$\displaystyle\,(L_{\rm ISCO}-L_{\rm BZ})\,\frac{1}{1+\eta_{\rm
J}}\dot{M}\,\mathbf{j}_{\rm BH}$ $\displaystyle-J_{\rm
BH}\bigg{\\{}\frac{\sin(\pi/7)}{\tau_{\rm GM}}(\mathbf{j}_{\rm
BH}\times\mathbf{j}_{\rm d})+\frac{\cos(\pi/7)}{\tau_{\rm
GM}}\big{[}\mathbf{j}_{\rm BH}\times(\mathbf{j}_{\rm BH}\times\mathbf{j}_{\rm
d})\big{]}\bigg{\\}}\,,$ (36)

and that of the disc is given by

$\displaystyle\dot{\mathbf{J}}_{\rm d}=$ $\displaystyle\,\dot{\mathbf{J}}_{\rm
d}^{\rm(acc)}-\dot{\mathbf{J}}_{\rm BH}^{\rm(BP)}\,,$ (37) $\displaystyle=$
$\displaystyle\,\dot{M}_{\rm in}\,\mathbf{L}_{\rm in}-L_{\rm
ISCO}\,\dot{M}\,\mathbf{j}_{\rm BH}$ $\displaystyle+J_{\rm
BH}\bigg{\\{}\frac{\sin(\pi/7)}{\tau_{\rm GM}}(\mathbf{j}_{\rm
BH}\times\mathbf{j}_{\rm d})+\frac{\cos(\pi/7)}{\tau_{\rm
GM}}\big{[}\mathbf{j}_{\rm BH}\times(\mathbf{j}_{\rm BH}\times\mathbf{j}_{\rm
d})\big{]}\bigg{\\}}\,.$ (38)

In terms of $\dot{M}$, the energy, mass and momentum fluxes into the jet are

$\dot{E}_{\rm J}=\frac{\epsilon_{\rm BZ}}{1+\eta_{\rm J}}\dot{M}c^{2}\,,$ (39)
$\dot{P}_{\rm J}=\frac{\sqrt{2\eta_{\rm J}\epsilon_{\rm BZ}}}{1+\eta_{\rm
J}}\,\dot{M}\,c\,,$ (40) $\dot{M}_{\rm J}=\frac{\eta_{\rm J}}{1+\eta_{\rm
J}}\dot{M}\,.$ (41)

##  3 Numerical implementation

Figure 1: A schematic overview of our sub-grid Blandford-Znajek jet model. On
the left we show how the sphere defined by the black hole smoothing length is
split into four sectors based on the direction of the $\alpha$-disc angular
momentum. The middle circle illustrates the jet injection in the context of
the wider simulation. The cylinder into which the jet is injected is outlined
in black and the white dashed line indicates the black hole smoothing length.
The circle on the right shows the angular momenta of the black hole and warped
$\alpha$-disc and how these correspond to the jet direction.

In the previous section we described the physical processes that govern our
sub-grid model and determined the equations that describe the evolution of the
system. In this section we detail how we have implemented this theoretical
model in our numerical framework.

###  3.1 AREPO

We have incorporated this new black hole accretion and feedback prescription
into the moving-mesh code arepo (Springel, 2010; Pakmor et al., 2016;
Weinberger et al., 2020). arepo solves the equations of hydrodynamics on a
quasi-Lagrangian, moving, unstructured mesh based on the Voronoi tessellation
of a set of points which move with the fluid with small corrections to ensure
regularisation of the cells. The code therefore inherits many benefits of
Lagrangian methods such as Galilean invariance and density dependent
resolution whilst also maintaining the high accuracy treatment of shocks and
instabilities that is characteristic of Eulerian codes. The gravitational
forces are calculated using a hierarchical oct-tree method (Springel, 2005;
Barnes & Hut, 1986) along with the PM method for long-range forces.

###  3.2 Gas inflow properties

At the beginning of each black hole timestep we first estimate the gas inflow
rate, $\dot{M}_{\rm in}$, and the specific angular momentum of this gas,
$\mathbf{L}_{\rm in}$, using a smooth particle hydrodynamic interpolation of
the local mass flux onto black hole, calculated within the black hole
smoothing length, $h_{\rm BH}$.

We do so by splitting the sphere defined by this smoothing length into four
sectors along the direction parallel to the angular momentum vector of the
$\alpha$-disc (as shown in Fig. 1). The reasoning behind splitting the sphere
into sectors will be explained shortly. The mass inflow rate is then
calculated for each sector individually according to

$\dot{M}_{{\rm
in},s}=-\int_{\mathcal{A}_{s}}\rho\mathbf{u}\cdot\mathrm{d}\mathbf{S}\approx-\mathcal{A}_{s}\mathcal{F}_{{\rm
BH},s}\,,$ (42)

where $s\in\\{1,2,3,4\\}$ indexes the sector, $\rho$ and $\mathbf{u}$ are the
density and velocity fields, respectively, and $\mathrm{d}\mathbf{S}$ is the
outward area element. $\mathcal{F}_{{\rm BH},s}$ is the numerical estimate of
the mass flux evaluated at the position of the black hole for this sector and
$\mathcal{A}_{s}$ is an effective area.

We approximate the mass flux in each sector by summing over all $N_{s}$ gas
cells that are radially inflowing and located within the relevant sector

$\mathcal{F}_{{\rm
BH},s}=\frac{\sum_{j=1}^{N_{s}}\rho_{j}u_{r,j}\,W(d_{j})}{\sum_{j=1}^{N_{s}}\,W({d_{j})}}\,,$
(43)

where $\rho_{j}$ is the mass density of the cell,
$u_{r,j}=(\mathbf{r}_{j}-\mathbf{r}_{\rm
BH})\cdot(\mathbf{u}_{j}-\mathbf{u}_{\rm
BH})/\left|\mathbf{r}_{j}-\mathbf{r}_{\rm BH}\right|$ is the radial velocity
of the cell and $W(d_{j})$ is a cubic spline kernel with compact support over
$h_{\rm BH}$ and $d_{j}=\left|\mathbf{r}_{j}-\mathbf{r}_{\rm BH}\right|/h_{\rm
BH}$.

We define the effective area to be $\mathcal{A}_{s}\equiv\pi r_{{\rm in},s}$
where $r_{{\rm in},s}$ is a characteristic inflow radius calculated for each
sector

$r_{{\rm
in},s}=\frac{\sum_{j=1}^{N_{s}}(3V_{j}/4\pi)^{1/3}\,W(d_{j})}{\sum_{j=1}^{N_{s}}\,W(d_{j})}\,,$
(44)

where $V_{j}$ is the volume of the gas cell.

Similarly, the specific angular momentum of the inflowing gas in each sector
is given by

$\mathbf{L}_{{\rm in},s}\equiv\frac{\dot{\mathbf{J}}_{{\rm
in},s}}{\dot{M}_{{\rm
in},s}}=\frac{\sum_{j=1}^{N_{s}}\rho_{j}u_{r,j}\,\mathbf{L}_{j}\,W(d_{j})}{\sum_{j=1}^{N_{s}}\rho_{j}u_{r,j}\,W({d_{j})}}\,,$
(45)

where $\mathbf{L}_{j}$ is the specific angular momentum of the cell.

We expect that inflowing gas will only be able to circularise and settle in
the disc if it has specific angular momentum that is smaller than that of the
outer edge of the disc

$L_{{\rm in},s}<L_{\rm d}(R_{\rm out})\approx\frac{J_{\rm d}}{M_{\rm d}}\,.$
(46)

We therefore check each sector (that contains inflowing material) individually
to ensure this criterion is met. If this condition is not satisfied then we
assume that there is no inflow in this sector and set $\dot{M}_{{\rm
in},s}=0$.

Carrying out this check for each sector individually allows us to better model
non-axisymmetric feeding of the black hole via coherent streams of material.
If accretion were to proceed in such a way, then calculating $\mathbf{L}_{\rm
in}$ using equation (45), but with the sum over all cells within $h_{\rm BH}$
may lead to the circularisation condition failing in scenarios where we would
expect the these streams to individually be able to reach the disc, due to the
nature of vector addition. We chose to use four sectors as we found that this
gave significant improvement in the accuracy of our inflow estimates whilst
still ensuring that each sector is sufficiently populated and that the
associated numerical overhead is minimised.

After this circularisation check has been carried out we are then able to
combine the values of $\dot{M}_{{\rm in},s}$ and $\mathbf{L}_{{\rm in},s}$ to
get a global values for the mass inflow rate, $\dot{M}_{\rm in}$, and the
specific angular momentum of inflowing material, $\mathbf{L}_{\rm in}$, by
summing over all sectors, with appropriate weightings. We store the values of
$\dot{M}_{{\rm in},s}$ for each sector to ensure that we drain the appropriate
amount of gas from each sector later in the timestep.

We now perform a series of further checks on $\dot{M}_{\rm in}$ and
$\mathbf{L}_{\rm in}$ to ensure that the assumptions of our model are
satisfied. We first verify that the specific angular momentum of the outer
edge of the $\alpha$-disc is larger than that required to maintain a circular
orbit at the ISCO (which holds true at all times in all simulations presented
in this work).

At this point we also cap $\dot{M}_{\rm in}$ to ensure that the predicted mass
of the accretion disc at the end of this timestep does not exceed the mass at
which we expect it to become self-gravitating (see equation (10)). Finally, we
limit $\dot{M}_{\rm in}$ to ensure that the disc remains sub-Eddington by
linearly interpolating equation (8) to ensure $f_{\rm Edd}<1$. If any of these
checks lead to $\dot{M}_{\rm in}$ being reduced then we also decrease the mass
inflow rate in each sector (where $\dot{M}_{{\rm in},s}\neq 0$) by the
relevant amount.

###  3.3 Black hole and disc mass evolution

The black hole and disc masses are then evolved by integrating equations (33)
and (34) over the timestep using a second-order predictor-corrector Heun’s
scheme. Letting the subscript ‘n’ indicate a value at the beginning of the
time-step; ‘temp’, a value after the predictor step and ‘n+1’, a value after
the corrector step, this integration procedure corresponds to

$\displaystyle M_{\rm BH,temp}=M_{\rm BH,n}+\bigg{(}\frac{1-\epsilon_{\rm
r}-\epsilon_{\rm BZ}}{1+\eta_{\rm J}}\bigg{)}\frac{f_{\rm Edd}}{\tau_{\rm
Salp}}\,M_{\rm BH,n}\Delta t\;,$ $\displaystyle M_{\rm d,temp}=M_{\rm
d,n}-\frac{f_{\rm Edd}}{\tau_{\rm Salp}}\,M_{\rm BH,n}\Delta t+\dot{M}_{\rm
in}\Delta t\;,$ $\displaystyle M_{\rm BH,n+1}=M_{\rm
BH,n}+\bigg{(}\frac{1-\epsilon_{\rm r}-\epsilon_{\rm BZ}}{1+\eta_{\rm
J}}\bigg{)}\frac{f_{\rm Edd}}{\tau_{\rm Salp}}\bigg{(}\frac{M_{\rm
BH,n}+M_{\rm BH,temp}}{2}\bigg{)}\Delta t\;,$ $\displaystyle M_{\rm
d,n+1}=M_{\rm d,n}-\frac{f_{\rm Edd}}{\tau_{\rm Salp}}\bigg{(}\frac{M_{\rm
BH,n}+M_{\rm BH,temp}}{2}\bigg{)}\Delta t+\dot{M}_{\rm in}\Delta t\;.$ (47)

###  3.4 Black hole and disc angular momentum evolution

The angular momenta of the black hole and disc evolve according to equations
(2.5) and (2.5). We implement this numerically by splitting the calculation
into two steps. In the first step we evolve the $\alpha$-disc due to inflow
from the surroundings and the loss of material that mass-loads the jet from
the ISCO. Simultaneously we evolve the black hole angular momentum according
to predictions of the Blandford-Znajek jet model

$\mathbf{J}_{\rm d,temp}=\mathbf{J}_{\rm d,n}+\mathbf{L}_{\rm in}\dot{M}_{\rm
in}\Delta t-\frac{\eta_{\rm J}}{1+\eta_{\rm J}}L_{\rm ISCO}\dot{M}\Delta
t\;\mathbf{j}_{\rm BH,n}\;,$ (48) $\mathbf{J}_{\rm BH,temp}=\mathbf{J}_{\rm
BH,n}-\frac{1}{1+\eta_{\rm J}}L_{\rm BZ}\;\dot{M}\;\Delta t\;\mathbf{j}_{\rm
BH,n}\,,$ (49)

where we are now using the subscript ‘temp’ to indicate the angular momenta
after this step. The factor of $\eta_{\rm J}/(1+\eta_{\rm J})$ in equation
(48) comes about due to the fact that we are only considering the material
that is drawn up into the jet from the inner edge of the disc, $\eta_{\rm
J}/(1+\eta_{\rm J})\;\dot{M}\;\Delta t$. Evolution due to the mass that flows
onto the black hole from the inner edge of the disc, $1/(1+\eta_{\rm
J})\;\dot{M}\;\Delta t$, is included in the subsequent step. We choose to
break down the calculations in this way because the total angular momentum of
the system after equations (48) and (49) have been evaluated, $\mathbf{J}_{\rm
cons}=\mathbf{J}_{\rm BH,temp}+\mathbf{J}_{\rm d,temp}$ is now conserved.

Having determined $\mathbf{J}_{\rm cons}$, we now begin the second step of the
calculation and evolve the disc and black hole due to the accretion of
material from the inner edge of the disc and due to mutual Bardeen-Petterson
torques.

The Bardeen-Petterson effect does not alter the magnitude of the black hole
angular momentum and also predicts that the inner disc will be
(counter-)aligned with the spin of the black hole. This means that the final
magnitude of the black hole angular momentum is affected by accretion alone

$J_{\rm BH,n+1}=J_{\rm BH,temp}+\frac{1}{1+\eta_{\rm J}}L_{\rm
ISCO}\dot{M}\Delta t\,.$ (50)

From this we calculate the predicted final spin of the black hole and cap it
to 0.998 if necessary.

Then we check to see if the black hole mass exceeds the warping mass (see
equation (16)). If not then the black hole angular momentum direction is
evolved using a predictor-corrector method according to

$\displaystyle\dot{\mathbf{j}}_{\rm BH,n+1}=$ $\displaystyle-
K_{1}(\mathbf{j}_{\rm BH,temp}\times\mathbf{j}_{\rm d,temp})$ $\displaystyle-
K_{2}\big{[}\mathbf{j}_{\rm BH,temp}\times(\mathbf{j}_{\rm
BH,temp}\times\mathbf{j}_{\rm d,temp})\big{]}\;,$ (51)

which applies the Bardeen-Petterson torque. If, however, the black hole mass
does exceed the warping mass then we assume that the angular momentum
directions of the black hole and disc align instantaneously with that of the
total angular momentum, as discussed in Section 2.2.

Finally, the disc angular momentum is updated to satisfy total angular
momentum conservation

$\mathbf{J}_{\rm d,n+1}=\mathbf{J}_{\rm cons}-\mathbf{J}_{\rm BH,n+1}\,.$ (52)

This final step incorporates the evolution of the disc due to the Bardeen-
Petterson effect and due to the disc material accreted by the black hole.

###  3.5 Mass draining

####  3.5.1 Sector based draining

Mass is then drained from each gas cell in the simulation that has been
determined to have contributed to the inflow onto the black hole-disc system
during this timestep, as described in Section 3.2. For a cell in a sector that
has inflow and is itself radially inflowing, we drain mass

$\Delta m_{i}={\rm max}\Bigg{(}\frac{m_{i}}{\sum_{j=0}^{N_{s}}m_{j}}\;\Delta
M_{{\rm in},s}\;,\;0.99\,m_{i}\Bigg{)}\,,$ (53)

where $i$ indexes the properties associated with an individual gas cell and
$\Delta M_{{\rm in},s}=\dot{M}_{{\rm in},s}\,\Delta t$.

####  3.5.2 Black hole dynamical mass evolution

In the simulation, the dynamical mass associated with the black hole particle
corresponds to that of the black hole-disc system as a whole, i.e the sum of
the black hole mass and the accretion disc mass. We update the dynamical mass
in accordance with this.

It should be highlighted here that black hole mass evolution is driven by
fluxes of binding energy associated with gas that falls across its horizon and
not the rest mass of this gas. This means that, in accordance with the
predictions of general relativity, black hole accretion should not explicitly
conserve mass, but rather energy. This lack of rest mass conservation is
particularly obvious for our sub-grid model in situations where there is no
net inflow of gas onto the system. In this scenario, no gas would be drained
from the simulation but the dynamical mass of the system may change due to the
accretion of material onto the black hole from the $\alpha$-disc.

###  3.6 Jet injection

We implement jet injection using a procedure similar to that of the ‘kinetic
energy conserving jet’ presented in Bourne & Sijacki (2017) which has been
adapted for our Blandford-Znajek model.

We inject the jet, on scales that can be resolved, into a cylinder of radius
$r_{\rm Jet}$ and height $2\,h_{\rm Jet}$ that is centred on the black hole.
The cylinder has axis that lies along the direction of the jet (i.e. the
direction of the black hole spin). $r_{\rm Jet}$ varies such that the total
mass within the jet cylinder, $M_{\rm cyl}$, is constant333We ensure a minimum
number of gas cells in each half of the jet cylinder meaning that the cylinder
mass could be larger than this target mass. In Section 3.7.2 we outline the
refinement techniques employed to ensure the jet cylinder is sufficiently
populated.. Throughout this work we choose our target jet cylinder mass to be
$M_{\rm cyl}=0.5\;{\rm M_{\odot}}$. $h_{\rm Jet}$ is determined by $r_{\rm
Jet}$ such that the aspect ratio of the cylinder is fixed, giving a constant
jet opening angle, $\theta_{\rm Jet}$, that satisfies

$\tan\bigg{(}\frac{\theta_{\rm Jet}}{2}\bigg{)}=\frac{r_{\rm Jet}}{h_{\rm
Jet}}\,.$ (54)

Following on from Bourne & Sijacki (2017), we choose $r_{\rm Jet}/h_{\rm
Jet}=3/2$ which means that the cylinder has total volume $(4\pi/3)\ r_{\rm
Jet}^{3}$. The cylinder is split into two halves, the north and the south
(each with height $h_{\rm Jet}$), and half of the required integrated energy,
mass and momentum flux is injected into gas cells in the ‘north’ and the other
half into those cells in the ‘south’ (see Fig. 1). The increase in mass of a
cell in the jet cylinder is given by

$\Delta m_{i}=\frac{1}{2}\dot{M}_{\rm J}\Delta t\;\frac{m_{i}\;W_{\rm
Jet}(r_{i},z_{i})}{M_{\rm weight}}\,,$ (55)

where $m_{i}$ is the initial mass of the cell, $\dot{M}_{\rm J}$ is as
determined by equation (41), the factor of a half accounts for the two halves
of the jet cylinder and $M_{\rm weight}$ is the weighted sum of the masses of
all gas cells in the relevant half-cylinder of the jet

$M_{\rm weight}=\sum^{\rm half-cyl}_{j}m_{j}\;W_{\rm Jet}(r_{j},z_{j})\,,$
(56)

where $r_{i}$ is the (cylindrical) distance of the gas cell from the jet axis
and $z_{i}$ is its height above the black hole equatorial plane. The kernel
function has the form

$W_{\rm Jet}(r,z)\propto\exp{\Bigg{(}-\frac{r^{2}}{2r_{\rm
Jet}^{2}}\Bigg{)}}\left|z\right|\,.$ (57)

Momentum injection is carried out by giving gas cells a kick along the jet
axis with magnitude

$\Delta p_{i}=\sqrt{2(m_{i}+\Delta m_{i})(E^{\rm kin}_{i}+\Delta E_{i}^{\rm
tot})}-\left|\mathbf{p}_{i}\right|\,,$ (58)

where $\mathbf{p}_{i}$ and $E^{\rm kin}_{i}$ are the initial momentum and
kinetic energy of the cell and $\Delta E^{\rm tot}_{i}$ is the desired change
in energy of the cell

$\Delta E^{\rm tot}_{i}=\frac{1}{2}\dot{E}_{\rm J}\,\Delta
t\;\frac{m_{i}\;W_{\rm Jet}(r_{i},z_{i})}{M_{\rm weight}}\,,$ (59)

where $\dot{E}_{\rm J}$ is as in equation (39).

This single momentum kick, however, does not necessarily ensure that the total
energy of the cell changes by $\Delta E_{i}^{\rm tot}$, as required by energy
conservation, unless the momentum vector of the cell was initially is parallel
to the jet axis. We therefore inject thermal energy, $\Delta\,E_{i}^{\rm
therm}$, to ensure energy conservation, where

$\Delta\,E_{i}^{\rm therm}=E^{\rm kin}_{i}+\Delta E_{i}^{\rm
tot}-\frac{1}{2}\frac{\left|\mathbf{p}_{i}+\Delta p_{i}\,\mathbf{j}_{\rm
BH}\right|^{2}}{(m_{i}+\Delta m_{i})}\,.$ (60)

It should be noted here that the Blandford-Znajek model predicts that the
black hole spins down as the jet is launched. In order to conserve angular
momentum we should also be giving the cells in the jet a kick in the toroidal
direction. We neglect this effect, since we expect it to correspond to
momentum kicks that are significantly smaller than those predicted by equation
(58).

###  3.7 Refinement

Our model makes considerable use of arepo’s ability to refine and de-refine
cells based on nearly arbitrary criteria which allows for high resolution in
specific regions of interest whilst maintaining lower resolution where not
needed. We employ refinement schemes that specifically target the region close
to the black hole, the jet injection cylinder and the jet lobes. Together
these allow us to accurately track the mass and angular momentum flows in the
vicinity of the black hole and follow the propagation of the jet at
sufficiently high resolution.

####  3.7.1 Black hole refinement

To obtain an accurate estimation the mass and angular momentum inflow rates,
we make significant use of the super-Lagrangian refinement scheme outlined in
Curtis & Sijacki (2015) which adaptively increases the mass and spatial
resolution in the region surrounding the black hole.

Within a chosen refinement radius, $R_{\rm ref}$, cells are forced to have
radii444When we refer to the cell radius, we mean the radius of a sphere with
volume equal to the cell volume: $R_{i}=\big{(}3V_{i}/4\pi\big{)}^{1/3}$.,
$R_{i}$, between two values, $R_{\rm min}(d)$ and $R_{\rm max}(d)$ which
increase linearly with distance from the black hole, $d$. Specifically, for a
cell with $d_{i}\leq R_{\rm ref}$, we ensure

$R_{\rm min}(d_{i})<R_{i}<R_{\rm max}(d_{i})\;,$ (61)

where

$\displaystyle R_{\rm min}(d)=\frac{d}{R_{\rm ref}}\frac{\big{(}R_{\rm
cell}^{\rm max}-R_{\rm cell}^{\rm min}\big{)}}{C}+\frac{R_{\rm cell}^{\rm
min}}{C}\;,$ (62) $\displaystyle R_{\rm max}(d)=\frac{d}{R_{\rm
ref}}\big{(}R_{\rm cell}^{\rm max}-R_{\rm cell}^{\rm min}\big{)}+R_{\rm
cell}^{\rm min}\;.$ (63)

Here, $C$, $R_{\rm cell}^{\rm max}$, and $R_{\rm cell}^{\rm min}$ are free
parameters to be specified. It should be noted that our choices for these
parameters differ from those suggested in Curtis & Sijacki (2015) as theirs
were chosen with the aim of resolving the Bondi radius. The simulations we
perform here are at a significantly higher spatial resolution and we are often
resolving the Bondi radius before this super-Lagrangian refinement is applied.

####  3.7.2 Jet refinement

Following the propagation of a jet at sufficiently high resolution presents
its own challenges and requires additional refinement schemes. We therefore
make use of further refinement prescriptions, outlined in Bourne & Sijacki
(2017). We briefly describe the relevant details of these schemes.

Firstly, from the jet axis out to a (cylindrical) radius of $\gamma r_{\rm
Jet}$, the cell masses are forced below a spatially dependent maximum mass
$m_{\rm cell}^{\rm max}(r)$ which ensures that the jet cylinder is
sufficiently well populated. $m_{\rm cell}^{\rm max}(r)$ decreases as one gets
closer to the jet axis according to

$m_{\rm cell}^{\rm max}(r)=\Bigg{[}(\alpha-\beta)\bigg{(}\frac{r}{r_{\rm
Jet}}\bigg{)}^{\kappa}+\beta\Bigg{]}M_{\rm cyl}\,,$ (64)

where r is the (cylindrical) distance to the jet axis, $\alpha$, $\beta$ and
$\gamma$ are parameters to be chosen and $\kappa$ is given by

$\kappa=\frac{1}{\ln(\gamma)}\ln\bigg{(}\frac{1-\beta}{\alpha-\beta}\bigg{)}\,,$
(65)

meaning that

$m_{\rm cell}^{\rm max}\big{(}\gamma\,r_{\rm Jet}\big{)}=M_{\rm cyl}\,.$ (66)

Both this refinement scheme and the black hole refinement scheme detailed
above only act within the $R_{\rm ref}$ (which always encloses the jet
cylinder) and so are localised to the vicinity of the black hole. Without
further refinement schemes, cells into which jet energy, momentum and mass are
injected are propelled out of this region and will then rapidly de-refine due
to arepo’s mass refinement criterion (Springel, 2010; Weinberger et al.,
2020).

We identify jet material using a passive tracer which tracks the jet material,
$m_{\rm i,J}$, in each cell. We set the tracer mass of cells in the jet
cylinder (into which jet material is injected) to the cell mass, $m_{i}$, and
then the tracer is advected with the gas. We make use of this ‘jet tracer’ to
ensure that there is sufficient resolution in the jet lobes by refining cells
with a jet-fraction ($f_{i,{\rm J}}\equiv m_{i,{\rm J}}/m_{i}$) greater than
$0.01$ if their volume satisfies

$V_{i}>\big{[}1-\log_{10}(f_{i,{\rm J}})\big{]}\,V_{\rm J}^{\rm max}\;,$ (67)

where $V_{\rm J}^{\rm max}$ is a parameter to be chosen.

We also employ further jet refinement schemes to stop cells de-refining if the
gradients between neighbouring cells are too steep (Bourne & Sijacki, 2017)
and to limit the volume of cells as well as the gradients between neighbouring
cells (Weinberger et al., 2017b).

###  3.8 Timesteps

The assumptions inherent in the Shakura & Sunyaev (1973) thin disc model and
those in our simplified model of the Bardeen-Petterson effect introduce
timescales that need to be resolved. We therefore ensure that the black hole
timestep satisfies

$\Delta t\leq 0.1\,{\rm min}\bigg{(}\tau_{\rm GM},\,\frac{M_{\rm
d}}{\dot{M}_{\rm d}}\bigg{)}\,,$ (68)

where $M_{\rm d}/\dot{M}_{\rm d}$ approximates the disc-draining time. Often,
however, further constraints on the black hole timestep that are not specific
to this sub-grid model force it to lower values (for further discussion see
Fiacconi et al. (2018)).

##  4 Simulation setup

Gas funnelled towards the centre of a galaxy due to secular evolution, large-
scale inflows or mergers may settle into a rotationally supported disc
surrounding the black hole. Indeed, such circumnuclear discs have been found
to be ubiquitous in interacting galaxies and those with signatures of recent
mergers (Medling et al., 2014; Imanishi & Nakanishi, 2013; Downes & Solomon,
1998). Observations have also shown that circumnuclear discs are
preferentially found in Seyferts compared to gas-rich spirals with no nuclear
activity (Hicks et al., 2013; García-Burillo & Combes, 2012; Davies et al.,
2007). These circumnuclear discs have radii $\sim 100$ pc and masses
$10^{8-9}\;{\rm M_{\odot}}$ (Chou et al., 2007; Schinnerer et al., 1999).
Motivated by the possibility of these discs residing at the centres of radio-
loud Seyferts we simulate an isolated circumnuclear disc embedded within a
stellar bulge and a single black hole at the centre. We surround the disc with
warm, diffuse gas, with properties analogous to that found in the CGM which
provides a medium into which the jet can propagate.

Modelling just the circumnuclear disc and neglecting any effects of the
surrounding galaxy, under the assumption that the circumnuclear disc is long
lived and in equilibrium, provides an ideal setup to test our model due to its
small scale nature as it allows us to focus our computational resources on the
jet propagation. Furthermore, this allows rapid evaluation of the model and
assessment of our parameter choices whilst still resolving the relevant
spatial and time scales.

###  4.1 Initial conditions

The initial placement of the mesh-generating points for the gas in the
circumnuclear disc and collisionless particles in the stellar bulge are
generated using a modified version of the procedure outlined in Springel et
al. (2005).

The stellar bulge is modelled using a Hernquist density profile (Hernquist,
1990)

$\rho_{*}(r)=\frac{M_{*}}{2\pi}\frac{r_{\rm s}}{r\,(r+r_{\rm s})^{3}}\,,$ (69)

where the bulge mass is $M_{*}=5\times 10^{8}\;{\rm M_{\odot}}$ and the scale
radius is $r_{\rm s}=100\;{\rm pc}$.

The gas disc is initialised with an exponential surface density profile

$\Sigma_{\rm gas}(R,z)=\frac{M_{\rm gas}}{2\pi\,h^{2}}\exp(-R/h)\,,$ (70)

where $M_{\rm gas}=10^{8}\;{\rm M_{\odot}}$ is the disc mass and $h\approx
70\;{\rm pc}$ is its scale-length. The initial vertical structure of the gas
disc is determined such that it is in hydrostatic equilibrium and is initially
isothermal with temperature $2\times 10^{4}\;{\rm K}$.

We sample the stellar bulge with $2\times 10^{6}$ collisionless particles of
mass $m_{*}=250\;{\rm M_{\odot}}$ and with gravitational softening length
$\epsilon_{*}=5\;{\rm pc}$. The gas disc is initialised with $4\times 10^{5}$
mesh generating points with target mass $m_{\rm t}=250\;{\rm M_{\odot}}$ and
gravitational softening length $\epsilon_{\rm gas}=0.5\;{\rm pc}$.

At the centre we place a black hole of mass $M_{\rm BH}=10^{6}\;{\rm
M_{\odot}}$ and gravitational softening length $\epsilon_{\rm BH}=5\;{\rm
pc}$. This black hole mass is consistent with those found in late type
galaxies with comparable stellar contents (Greene et al., 2019).

We then add additional mesh generating points to fill the box with a diffuse,
hot CGM. The gas cells are placed at rest, with temperature $10^{7}$ K and
with their positions sampled from a uniform distribution to give a constant
density. We consider the case where this CGM material has density $\rho_{\rm
CGM}=10^{-27}\;{\rm g\,cm^{-3}}$ (chosen such that the circumnuclear is in
approximate pressure equilibrium with the CGM) which we refer to as the
‘standard’ CGM and we also consider the case where the CGM is two orders of
magnitude more dense, i.e. $\rho_{\rm CGM}=10^{-25}\;{\rm g\,cm^{-3}}$, which
we refer to as the ‘dense’ CGM.

These initial conditions are then evolved for $200\;{\rm Myrs}$ (approximately
$4$ orbital periods of the disc) to allow transient features in the disc to
dissipate and for the disc to come into pressure equilibrium with the CGM. The
increased pressure in the ‘dense’ CGM causes noticeable flattening of the
circumnuclear disc and more instabilities are seen along the disc-CGM
interface, however the disc is still able to reach quasi-steady state during
the relaxation period.

###  4.2 Additional tracers

In these simulations we would like to be able to track how much disc material
is physically entrained by the jet. When generating the initial conditions, we
therefore initialise a passive Lagrangian tracer, $m_{i,{\rm disc}}$, in all
gas cells in the circumnuclear disc, which we refer to as the ‘disc tracer’.
The value of this disc tracer is then set to zero in all cells into which jet
mass, momentum and energy are injected meaning that it tracks pure disc
material.

Since our disc and jet tracers are independent, we obtain an effective CGM
tracer for each cell via

$m_{i,{\rm CGM}}=m_{\rm i}-m_{i,{\rm J}}-m_{i,{\rm disc}}\,.$ (71)

We use this tracer to asses the levels of mixing and entrainment of CGM
material in the jet lobes.

###  4.3 Additional refinement

As discussed in Section 3.7, the mass refinement scheme in arepo maintains gas
cells at a chosen target mass, which we choose to be $m_{\rm t}=250\;{\rm
M_{\odot}}$ as this allows us to maintain sufficiently high resolution in the
circumnuclear disc. Since the velocity of propagation of jet material is often
supersonic, cells in the CGM are frequently unable to refine before the jet
material reaches them which is particularly problematic if we maintain
relatively coarse spatial resolution in the CGM. This can lead to increased
numerical mixing at the jet-CGM interface as well as instabilities being
washed out and the jet morphology and propagation being altered significantly
(see e.g. Bourne & Sijacki, 2017).

Due to the density contrast between the circumnuclear disc and the CGM it is
not necessarily the case that a global target mass $m_{\rm t}=250\;{\rm
M_{\odot}}$ will yield sufficiently high resolution in the CGM to mitigate
these problems. In the ‘dense’ CGM we find that we are, in fact, able to
achieve sufficient resolution in the CGM with just a global target mass of
$250\;{\rm M_{\odot}}$. In the ‘standard’ CGM, however, this is not the case
and we, therefore, implement an additional refinement scheme in any simulation
that has a ‘standard’ density CGM.

Since it is numerically infeasible to increase the resolution in the entire
CGM, we enforce a spatially dependent, piecewise-defined target mass that is
specifically applied to the gas cells in the CGM555Note that we do not apply
this refinement scheme to cells in the jet. If the refinement scheme were only
applied to cells in the disc then any such disc cell that enters the jet
cylinder and becomes jet material will suddenly have its target mass decreased
and we find that the resulting rapid refinement lead to artificial features in
the jet., targeting the increased resolution to a central cylindrical region
(with the cylinder axis parallel to the $z$-axis) in which the majority of the
jet-CGM interaction is expected to occur. Specifically, we define the target
mass for CGM cells at a cylindrical radius, $R$, to be

$m_{\rm t}^{\rm CGM}(R)=\begin{cases}f_{\rm ref}\,m_{\rm t}&R<R^{\rm CGM}_{\rm
ref}\\\ (a_{\rm ref}R+b_{\rm ref})\,m_{\rm t}&R^{\rm CGM}_{\rm ref}<R<N_{\rm
ref}\,R^{\rm CGM}_{\rm ref}\\\ m_{\rm t}&R>N_{\rm ref}\,R^{\rm CGM}_{\rm
ref}\end{cases}$ (72)

where $f_{\rm ref}$, $N_{\rm ref}$ and $R^{\rm CGM}_{\rm ref}$ are free
parameters.
The parameters $a_{\rm ref}$ and $b_{\rm ref}$ are chosen such that the
$m_{\rm t}^{\rm CGM}(R)$ is continuous in $R$, i.e.

$\displaystyle a_{\rm ref}=\frac{1}{R_{\rm ref}}\bigg{(}\frac{1-f_{\rm
ref}}{N_{\rm ref}-1}\bigg{)}\,,$ (73) $\displaystyle b_{\rm
ref}=\bigg{(}\frac{f_{\rm ref}\;N_{\rm ref}-1}{N_{\rm ref}-1}\bigg{)}\,.$ (74)

##  5 Jets with fixed power and direction: Results

When implementing the sub-grid black hole accretion and Blandford-Znajek jet
model, detailed in the previous sections, we first need to validate it against
analytic expectations to verify the physical nature of our jets. Due to the
complexity of the full sub-grid model, however, we first carry out a series of
high-resolution studies in which the jet power and direction are fixed. Taking
away these degrees of freedom makes for a more straightforward comparison with
analytical models and creates a benchmark against which we can compare the
results of runs that use our full sub-grid model. This simplified setup also
allow us to more cleanly isolate how jets are expected to interact with the
surrounding circumnuclear disc and the CGM.

In these idealised simulations we fix the jet power to a fraction of the
Eddington luminosity of the black hole, $L_{\rm Edd}$666Here $L_{\rm
Edd}\approx 1.26\times 10^{44}\;{\rm erg\,s^{-1}}$ refers to the Eddington
luminosity of the initial black hole mass, $M_{\rm BH}=10^{6}\;{\rm
M_{\odot}}$, and disregards the fact that the Eddington luminosity of the
black holes will change very slightly throughout the course of the simulation
as the black hole accretes., and fix the direction of the jet to lie parallel
to the z-axis.

Our simulation suite consists of six runs in total. We consider jets with
three different powers: $10^{-4}L_{\rm Edd}$ \- the ‘low’ power run,
$10^{-3}L_{\rm Edd}$ \- the ‘medium’ power run, and $10^{-2}L_{\rm Edd}$ \-
the ‘high’ power run. These jets are then launched from the centre of the disc
and we consider both the ‘standard’ and ‘dense’ CGM. We also fix the sub-grid
velocity to $v_{\rm J}=10^{5}\;{\rm km\,s^{-1}}$. Doing so removes all the
spin-dependence of the jet energy, momentum and mass fluxes; these are
constant and fully specified by equations (28) and (29). A full list of the
parameters of these six simulations can be found in Table 1.

These runs are of a similar nature to those in Bourne & Sijacki (2017) in that
they utilise the ‘kinetic’ jet model and the jets have fixed direction. In our
simulations, however, the jet power is fixed to some fraction of the initial
Eddington luminosity of the black hole whereas the jet power in Bourne &
Sijacki (2017) may change slightly as the black hole grows (note however that
in those simulations the black hole growth was insignificant). In addition, we
are probing significantly smaller scales777In Bourne & Sijacki (2017), the
target jet cylinder mass is $10^{4}\;{\rm M_{\odot}}$ whereas in this work it
is $0.5\;{\rm M_{\odot}}$. and, crucially, we are launching jets from a black
hole embedded in a circumnuclear disc, whereas Bourne & Sijacki (2017)
launched jets from black holes embedded in a hydrostatic atmosphere.

Table 1: Properties of the simulations in which the jet power and direction are fixed. The first and second columns show textual identifiers that we use throughout the results section when we refer to a specific simulation (e.g. the ‘high’ power jet launched into the ‘dense’ CGM). The power of the jet is shown in the third column. The fourth column corresponds to the density of the CGM. The fifth and sixth columns show the final time of each simulation and the length of the northern lobe of the jet at this time, respectively. Power label | CGM label | Jet power | CGM density | Final time | Jet length at final time
---|---|---|---|---|---
|  | [${\rm ergs\,s^{-1}}$ : $L_{\rm Edd}$] | [${\rm g\,cm^{-3}}$] | [Myrs] | [pc]
‘high’ | ‘standard’ |  $1.26\times 10^{42}$ : $10^{-2}$ | $10^{-27}$ | $0.43$ | $2964$
‘medium’ | ‘standard’ |  $1.26\times 10^{41}$ : $10^{-3}$ | $10^{-27}$ | $0.97$ | $2989$
‘low’ | ‘standard’ |  $1.26\times 10^{40}$ : $10^{-4}$ | $10^{-27}$ | $2.54$ | $2995$
‘high’ | ‘dense’ |  $1.26\times 10^{42}$ : $10^{-2}$ | $10^{-25}$ | $1.89$ | $2996$
‘medium’ | ‘dense’ |  $1.26\times 10^{41}$ : $10^{-3}$ | $10^{-25}$ | $5.87$ | $2999$
‘low’ | ‘dense’ |  $1.26\times 10^{40}$ : $10^{-4}$ | $10^{-25}$ | $10.00$ | $2120$
Figure 2: From top to bottom, each row shows slices through the computational
domain for the ‘high’, ‘medium’ and ‘low’ power jet runs in the ‘standard’ CGM
case. The first and second columns show the temperature and density fields,
respectively. The third column shows a pressure diagnostic (which is described
in Section 5.1). The fourth column shows the thermal pressure field and then
the fifth column shows the $y$/$z$ component of the vorticity field in the
top/bottom of each slice. Finally, in the sixth column, the slices indicate
the internal Mach number of the gas, overlaid with streamlines of the velocity
field in the $x$-$z$ plane, where the width of the streamline scales with the
logarithm of the absolute velocity field. Each slice was created at the time
at which the length of the northern lobe of the jet first exceeds $2000\;{\rm
pc}$. This instant in time is then indicated in the top right-hand corner of
the relevant temperature slice. Figure 3: From top to bottom, each row shows
slices through the computational domain for the ‘high’, ‘medium’ and ‘low’
power jet runs in the ‘dense’ CGM case. The first and second columns show the
temperature and density fields, respectively. The third column shows a
pressure diagnostic (which is described in Section 5.1). The fourth column
shows the thermal pressure field and then the fifth column shows the $y$/$z$
component of the vorticity field in the top/bottom of each slice. Finally, in
the sixth column, the slices indicate the internal Mach number of the gas,
overlaid with streamlines of the velocity field in the $x$-$z$ plane, where
the width of the streamline scales with the logarithm of the absolute velocity
field. Each slice was created at the time at which the length of the northern
lobe of the jet first exceeds $2000\;{\rm pc}$. This instant in time is then
indicated in the top right-hand corner of the relevant temperature slice.

###  5.1 Basic jet properties

We begin this section by providing a qualitative overview of the main
morphological properties of these fixed power jets and discuss the physical
processes that are driving the formation of these features. We perform this
analysis with the aid of Figs. 2 and 3 which show slices through the
computational domain of various diagnostic quantities for the runs in the
‘standard’ and ‘dense’ CGM, respectively. From top to bottom, each row shows
the relevant slices for the ‘high’, ‘medium’ and ‘low’ power jet runs. Since
the propagation timescales of the jets varies significantly from run-to-run
and we wish to make comparisons between the jets on the same spatial scale, we
select the snapshot at which the length of the northern lobe of each jet first
exceeds $2000\;{\rm pc}$ and generate these slices at this time. The relevant
time for each run is then indicated in the top right-hand corner of the
corresponding temperature slice.

In each figure, the first and second columns show the temperature and density
fields, respectively. To better identify the regions in which ram pressure or
thermal pressure dominate, we define a pressure ratio

$P_{\rm ratio}=\frac{P_{\rm ram}-P_{\rm therm}}{P_{\rm ram}+P_{\rm therm}}\,,$
(75)

where $P_{\rm ram}=\rho u^{2}$ is the ram pressure and $P_{\rm therm}$ is the
thermal pressure (Bourne & Sijacki, 2017). Slices of this quantity are shown
in the third column. In the fourth column we show the thermal pressure field
and then in the fifth, we show the $y$/$z$ component of the vorticity field in
the top/bottom of each slice, where the red and blue colour maps indicate
regions where the relevant vorticity component has a positive or negative
sign, respectively. Finally, in the sixth column, the slices indicate the
internal Mach number of the gas, onto which we have overlaid streamlines of
the velocity field in the $x$-$z$ plane, where the width of the streamline
scales with the logarithm of the absolute velocity field.

Unsurprisingly, in all runs we can clearly identify the centrifugally
supported circumnuclear disc from which the jet has been launched. In the
‘dense’ runs, however, the disc has a noticeably smaller scale-height and
scale-length as the higher pressure in this CGM has somewhat compressed the
disc gas. Both in the ‘standard’ and in the ‘dense’ case, when the jets first
become active they immediately encounter resistance from this cold, dense disc
material and drive strong shocks into the disc. In general, the ‘high’ power
jets are more effective at pushing aside the disc material and emerge into the
CGM with highly supersonic bulk velocities. The ‘low’ power jets, however, are
not as effective at overcoming the inertia of the disc material and, as a
result, are colder and slower, albeit still supersonic. This can be quantified
by considering the maximum Mach numbers within the jet material during the
first $0.1\;{\rm Myrs}$ after the jet is launched. For the high power jet we
find Mach numbers of $\sim 270$ and $\sim 160$ for the ‘standard’ and ‘dense’
jets, respectively, whereas those found in the ‘low’ power jets are
significantly lower: $\sim 50$ for the ‘standard’ run and $\sim 40$ for the
‘dense’ run, although these still indicate the presence of strong shocks. As
the jets leave the disc, some of the disc material is entrained in the
outflow, clearly identifiable in all runs by the dense regions surrounding the
jet channels and, indeed, close to the base of the jet this disc component
dominates the mass of the jets (see Sections 5.3 and 5.4 below for further
discussion).

In all runs, the launching of the jets also drives lateral shocks into the
disc which advance in an almost hourglass shape in the radial direction,
significantly altering the thermodynamic profiles of the inner disc, with the
higher power jets able to influence disc material out to larger radii. These
shocks rapidly weaken as they propagate outwards, eventually broadening into
sound waves, leaving the outer disc largely unperturbed for the duration of
these simulations. The interaction of these shocks and their reflections sets
up a complex pattern of interference which mediates the properties of the gas
in the vicinity of the black hole and as we do not allow the disc gas to
radiatively cool these waves are not able to damp efficiently.

Outside of the disc, the launching of the jet drives a strong horizontal shock
into the CGM with lateral velocity significantly higher than the velocity of
the shock travelling within the disc material. This leads to an almost peanut-
shaped shock front that advances into the CGM, particularly apparent in the
pressure maps. In all runs this shock front initially envelops the whole jet,
however, with time the strength of its horizontal advance diminishes and
eventually becomes subsonic, leaving only the head of the jet driving shocks
into the CGM. Over time these forward shocks also broaden into sound waves and
detach from the head of the jet. At the times at which the slices in Figs. 2
and 3 were made, the ‘medium’ and ‘low’ power jets in the ‘dense’ CGM are no
longer driving shocks into the CGM, whereas the ‘high’ power jet in the
‘dense’ CGM and all jets in the ‘standard’ CGM are still mostly enveloped by
shock fronts.

Turning now to the jets themselves, it is evident from inspection of Figs. 2
and 3 that the morphology and propagation of these jets are sensitive both to
the power of the jet and to the density of the surrounding medium. These
figures also clearly show that, on the time and length scales considered here,
the jets in the ‘dense’ CGM are particularly sensitive to changes in jet power
as they display a much more diverse range of morphological characteristics
than those launched into the ‘standard’ CGM which look comparatively similar.

In Fig. 3 the pressure and temperature slices highlight the series of
recollimation shocks that have formed along the jet axis in all the ‘dense’
CGM runs. These shocks form because, as the jets exit the disc, they are
initially overpressured with respect to the external CGM and expand laterally
until eventually they become under-pressured and are focused by the
surrounding CGM. Since perturbations at the jet channel boundary are
communicated by acoustic waves, the jet channel boundary oscillates as it
repeatedly overshoots the equilibrium configuration, forming the recollimation
shocks. Due to the lower CGM pressure in the ‘standard’ CGM runs the length-
scales over which these jets are recollimated are larger than the scales
considered in Fig. 2. We do, however, observe that the jet beams narrow
towards terminal shock as the recollimation process begins. At later times,
when these jets are eventually recollimated, they however do not display the
repeated pattern of shocks that is so obvious in the ‘dense’ CGM runs.

In the ‘dense’ CGM runs the jet beam gas alternately heats and is compressed
as it passes through the recollimation shocks and cools and expands as it
passes through the rarefaction fans. Once the first recollimation shock has
formed in each run, it remains approximately stationary with a spatial offset
from the black hole that increases with jet power, ranging from $\sim
150\;{\rm pc}$ in the ‘low’ power run to $\sim 400\;{\rm pc}$ in the ‘high’
power run. The fact that the jets in the ‘standard’ CGM are recollimated on
significantly larger length-scales is indicative that, in this setup, it is
the CGM density rather than the power of the jet that has a more significant
influence on the timescale of the recollimation process.

Upstream of this first recollimation shock, we find fluid instabilities
developing in the material surrounding the jet channel. However, it is
downstream where these instabilities come to dominate the flow structure where
they act to disrupt the jet beam and divert shocked jet material off-axis,
enhancing the levels of turbulence and mixing in the surrounding cocoon.
Indeed, as evidenced in Fig. 3, in the ‘low’ power jet these instabilities
have completely disrupted the jet beam whereas the ‘high’ power jet still
retains a ram pressure-dominated beam which terminates in a weaker shock,
close to the location of the forward shock being driven into the CGM.

As noted above, the morphology of the jets in the ‘standard’ CGM runs is
fairly similar, but there are some differences in the thermodynamic properties
that are worth noting. Looking specifically at the jet channel region, as we
move from low to high jet powers the temperature increases and the density
decreases in general. The decrease in temperature with jet power can be used
to understand why in Fig. 2 we find the highest internal Mach numbers in the
‘low’ power jet. Since this quantity corresponds to the ratio of the absolute
velocity of each gas cell to its sound speed, as the jet power decreases, the
lower sound speeds associated with lower temperatures are sufficient to offset
the decreasing velocities.

The temperature, density and pressure maps of the ‘standard’ CGM simulations
show evidence for a series of weak internal shocks that are propagating along
the jet channel as shells of material travelling with different velocities
collide. These velocity differences may be due to changes in the thermodynamic
properties of the injection region which lead to fluctuations in the initial
mass loading but we also expect that interactions with the circumnuclear disc
play a significant role here. While we do find evidence for similar internal
shocks in the jets in the ‘dense’ CGM runs these are harder to see in Fig. 3
due to the smaller width of the jet channels.

The narrowing of the jet channel towards the head of the jets in the
‘standard’ CGM runs leads to the formation of oblique shocks along the channel
boundary and the jets then terminate in a strong annular shock at the pole of
the jet axis in which the majority of the ram pressure dominated jet beam
material undergoes thermalisation. This then leads to the formation of a ‘hot
spot’ at the stagnation point of the flow as the beam material is essentially
brought to a halt relative to the downstream contact discontinuity. This
hotspot, particularly visible in the pressure and temperature maps, is more
obvious in the case of the ‘high’ power jet. We do, however, see evidence for
similar features in the ‘medium’ and ‘low’ power jets. The ram pressures of
the beam and shocked CGM felt by jet material passing through the terminal
shock causes it to be deflected laterally to regions of lower total pressure
and it forms a wide fan of backflowing material. This material, along with
that which has passed through the oblique shocks, feeds the low-density cocoon
of shocked jet material. Conversely, in the ‘dense’ CGM runs, the strongest
on-axis shock is the first recollimation shock through which jet material
passes and so the cocoon is fed due to the action of instabilities that have
been excited by the recollimation shocks which act in addition to the terminal
shock, where it exists.

For all three jet powers in the ‘dense’ CGM case there is a significant
enhancement in both components of the vorticity, particularly in the turbulent
cocoon but also in the shocked CGM enveloped by the bow shock. There is also
evidence for enhanced levels of vorticity in all of the ‘standard’ runs, with
greater levels of vorticity generation found for higher power jets. Vorticity
enhancements are expected to come about due to turbulence and g-modes. The
fact that we do not see evidence in any of the runs for significant vorticity
enhancements outside of the jet cocoon indicates that these jets are not able
to drive significant turbulence in the CGM (Bourne & Sijacki, 2017; Weinberger
et al., 2017b; Yang & Reynolds, 2016a; Reynolds et al., 2015) and that any
g-modes that have been excited are trapped within the cocoon for the duration
of our simulations (Reynolds et al., 2015).

###  5.2 Comparison with analytic predictions

Figure 4: Individual points correspond to a value measured from the simulation
data and the corresponding lines are the best-fit analytic model for the jet
length (top row), the velocity of the jet head (middle row) and $r_{\rm
ws}/\chi^{1/2}$ (bottom row). These are plotted as a function of $t-t_{0}$
where $t_{0}$ is the time at which the relevant jet has length: $Z_{\rm
J,0}=500\;{\rm pc}$. The ‘standard’ CGM runs are shown in the left column and
the ‘dense’ CGM runs are shown in the right column.

Having explored some of the salient features of these jets we now undertake a
more quantitative analysis and compare their evolution to analytic
predictions. We do not perform these comparisons with the aim of showing that
our jets are well described by a specific model since we expect that the
initial interaction with the circumnuclear disc will affect the long-term
evolution of the jet (at least to some extent) and, as highlighted in Section
5.1, significant instabilities can develop in certain cases. We perform these
comparisons, rather, with the aim of identifying regimes in which assumptions
inherent to the analytic model either hold or break down which can be used to
aid our understanding of the dominant processes affecting jet evolution.

Several works have developed analytic models that predict the time evolution
of the jet length (Harrison et al., 2018; Bromberg et al., 2011; Perucho &
Martí, 2007; Kaiser & Alexander, 1997; Falle, 1991; Begelman & Cioffi, 1989).
Here, we choose to provide a comparison against the seminal analytic model
detailed in Begelman & Cioffi (1989) which predicts the time evolution of the
length of the jet by balancing the thrust of the jet, $\Pi_{\rm J}$, with the
ram pressure force of a uniform density CGM888Begelman & Cioffi (1989)
actually describe the medium into which the jet is launched as the ICM. Since
the only relevant assumptions are that the medium has uniform density, we here
refer to it as CGM material, to remain consistent with our setup.:

$\rho_{\rm CGM}v_{\rm h}^{2}A_{\rm h}=\Pi_{\rm J}\,,$ (76)

where $A_{\rm h}$ is the cross-sectional area of the bow shock at the end of
the cocoon and $v_{\rm h}$ is the velocity of the head of the jet. Begelman &
Cioffi (1989) assume that the thrust of the jet is given by $\Pi_{\rm J}\sim
L_{\rm J}/v_{\rm J}$ where $L_{\rm J}$ is the jet power and $v_{\rm J}$ is the
jet velocity.

Bourne & Sijacki (2017) used this model to validate the jets in their
simulations, under the assumption that the jet thrust is given by $\Pi_{\rm
J}=\chi\dot{M}_{\rm J}v_{\rm J}/2$ where $\dot{M}_{\rm J}$ is the sub-grid jet
mass flux, $v_{\rm J}$ is now the sub-grid jet velocity and the factor of
$1/2$ comes about due to the momentum flux being split between the two halves
of the jet cylinder. The factor $\chi$ is a constant which takes into account
any momentum boost that occurs in the initial launching of the jet. Bourne &
Sijacki (2017) also expressed the cross-sectional area of the bow shock in
terms of a working surface radius: $A_{\rm h}=\pi r_{\rm ws}^{2}$ and, under
the above assumptions, found that their jets were a fairly good fit to the
analytic model if the working surface radius scales linearly with time.

It should be noted here that implicit in the choice that $\chi$ be a constant
is the assumption that the entire momentum flux at the base of the jet,
$\chi\dot{M}_{\rm J}v_{\rm J}/2$, is what is felt at the head of the jet. In
our analysis, we relax the assumption that $\chi$ is constant and instead of
assuming that $r_{\rm ws}$ scales linearly with time, we make a more general
assumption that $r_{\rm ws}/\chi^{1/2}=at+b$ for some constants $a$ and $b$.
This is consistent with the fitting procedure in Bourne & Sijacki (2017) but
alters the interpretation of the result. In our case, it should be made clear
that $r_{\rm ws}/\chi^{1/2}$ does not give us any information about the
behaviour of the working surface radius or the momentum flux at the head of
the jet as the two are degenerate.

To determine the analytic jet length evolution we obtain the velocity of the
jet head from equation (76)

$v_{\rm h}\approx\bigg{(}\frac{\chi\dot{M}_{\rm J}v_{\rm J}}{2\pi\rho_{\rm
CGM}}\bigg{)}^{1/2}\frac{1}{r_{\rm ws}}\,.$ (77)

Inserting $r_{\rm ws}/\chi^{1/2}=at+b$ then gives

$v_{\rm h}(t)\approx\bigg{(}\frac{\dot{M}_{\rm J}v_{\rm J}}{2\pi\rho_{\rm
CGM}}\bigg{)}^{1/2}\frac{1}{at+b}\,,$ (78)

and integrating this then gives the jet length

$Z_{\rm J}(t)\approx\bigg{(}\frac{\dot{M}_{\rm J}v_{\rm J}}{2\pi\rho_{\rm
CGM}}\bigg{)}^{1/2}\frac{1}{a}\ln\bigg{(}\frac{at+b}{at_{0}+b}\bigg{)}+Z_{\rm
J,0}\,,$ (79)

where $t_{0}$ is the time at which the jet length is $Z_{\rm J,0}$.

Since initial interactions with the circumnuclear disc (and differences in the
initial mass-loading between the north and south halves of the jet cylinder)
can lead to slight differences in the evolution of the north and south lobes
of a given jet (which are typically more pronounced at early times), we
restrict our analysis to the northern jet lobe. Further to this, since we know
this model is not valid for the initial evolution of the jets before they are
able to break out of the circumnuclear disc, for each jet we begin the
comparison at the time when its length has reached $Z_{\rm J,0}=500\;{\rm
pc}$. We then follow the evolution of each jet for a period of time that is
sufficient to be able to make robust comparisons with the analytic model
(specifically, until its length reaches $3000\;{\rm pc}$ or $2\;{\rm Myrs}$
has elapsed). After the length of the northern jet lobe has been calculated,
the velocity of the jet head is found by differencing the jet length between
snapshots and we then fit for the constants $a$ and $b$.

Figure 5: The positive flux of $z$-momentum in jet material as a function of
distance along the jet axis. The top and bottom rows show results from runs in
the ‘standard’ and ‘dense’ CGM cases, respectively. Each column, from left to
right, corresponds to the ‘high’, ‘medium’ and ‘low’ power jets. Individual
lines show the momentum flux profile for the relevant jet at a specific time
where the times have been chosen to sample the whole range for which the
analytic fitting procedure was carried out (see Fig. 4). Lighter lines
correspond to a larger value of $t-t_{0}$. For each line, the arrow on the
$x$-axis with the corresponding weight indicates the length of the jet.

The values measured from the simulation and the best fit analytic models are
displayed in Fig. 4. The top row shows the jet length, the second row shows
the velocity of the jet head and the third row shows the resulting value for
the quantity $r_{\rm ws}/\chi^{1/2}$. Each individual point corresponds a
value measured from the simulation data and the solid line shows the best-fit
analytic prediction. The results for the jets in the ‘standard’ CGM case are
shown in the left column and the jets in the ‘dense’ CGM case are shown on the
right999Whilst the fits were performed for each jet using the data from all
snapshots in the relevant time interval, for ease of viewing we have plotted
every third data point for the jet lengths..

Inspecting this figure it is immediately clear that, on the whole, the
analytic model provides a relatively good fit to the jets in the ‘standard’
CGM. These runs also display the expected scaling of the jet length and
velocity with jet power, with the largest jet head velocities found in the
‘high’ power jet. The jets in the ‘dense’ CGM, however, do not seem to be well
described by this model and show significantly more variability in the
measured velocity. Moreover, the ‘low’ power jet in the ‘dense’ CGM shows two
distinct phases in its evolution with the initial jet propagation stalling
around $1.17\;{\rm Myrs}$ after the jet turns on (which corresponds to
$t-t_{0}=0.49\;{\rm Myrs}$) after which its length stagnates. This behaviour
is obviously not well fit by just a single pair of values, ($a$, $b$), and so
the analytic fit we provide in Fig. 4 corresponds to a piece-wise function
which better reflects this bimodal behaviour. Whilst there is not such a clear
transition in the ‘medium’ and ‘high’ power jet runs, there are still
differences between their early and late time behaviour that have not been
well captured by the analytic model.

We further examine this by considering the positive flux of z-momentum in jet
material101010Here we define jet material to be all cells with $f_{\rm
J}>10^{-3}$ (see Section 5.3 for further discussion). Whilst this definition
will lead to cocoon material that is outflowing being included in the
calculation, we do not expect this to change our interpretation of the
momentum fluxes. as a function of distance along the jet axis, shown in Fig.
5. The first and second rows show the results from runs in the ‘standard’ and
‘dense’ CGM cases, respectively, with each column corresponding to the ‘high’,
‘medium’ and ‘low’ power jets from left to right. Each line then corresponds
to the momentum flux profile for the relevant jet at a specific time where the
times have been chosen to sample the whole range for which the analytic
fitting procedure was carried out, with lighter lines corresponding to a
larger values of $t-t_{0}$.

We have made use of arepo’s shock finder (Schaal & Springel, 2015) to
characterise the shocks within the jet material. On each line, the marker ‘x’
shows the location of the strongest shock in the outflowing jet material at
this time and the colour of the marker indicates the Mach number of this
shock. As expected, following the discussions in Section 5.1, the strongest
shock in the ‘standard’ CGM runs is the reverse shock at the head of the jet.
In the ‘dense’ CGM runs, before the recollimation shock forms, the strongest
shock is also the reverse shock at the head, but once the recollimation shock
forms, this first shock through which jet material passes is strongest.
Comparing the strength of the shocks we see that those in the ‘dense’ CGM runs
typically have lower Mach numbers than those in the ‘standard’ CGM.

For the ‘high’ and ‘medium’ power runs we see that, in general, the momentum
flux decreases along the length of the jet, but this decrease is much more
significant in the ‘dense’ CGM runs and it primarily occurs beyond the
recollimation shock, as material is thermalised and diverted off-axis.

We also see that for most jets that there is a peak in momentum flux at early
times which propagates towards the head before settling down somewhat. This
behaviour is akin to a transient in the initial conditions and occurs
primarily due to the gas in the vicinity of the black hole readjusting to the
launching of the jet. Specifically, as the density in the centre drops, the
action of the black hole refinement allows the jet cylinder mass to fall down
to its target mass which, in turn, increases the velocity of individual cells
kicked out of the jet cylinder. The communication of these higher velocities
to the head of the jet is what leads to the travelling peak seen in the
momentum flux profiles. At higher jet powers, regulation of the gas properties
in the vicinity of the black hole occurs on shorter timescales and so this
feature persists for longer in the ‘low’ power jets.

Crucially however, in all runs but particularly those in the ‘dense’ CGM, the
momentum flux at the head of the jet is not constant in time; rather it seems
to decrease. It is therefore likely that the break down of the analytic model
can be attributed to the assumption of jet momentum flux conservation and, in
particular, in the ‘dense’ CGM case this occurs due to the presence of strong
recollimation shocks which thermalise jet material and the ensuing onset of
significant instabilities which divert jet beam material off-axis.
Furthermore, the fact that the momentum flux at the head of the jet is time-
variable justifies our interpretation that $\chi$ should be characterised as a
time dependent quantity that corresponds to the momentum flux felt by the head
of the jet in units of the sub-grid momentum flux. In the case where multiple
shocks and turbulence are acting, this is a complex non-linear function that
is not well described by a simple analytical model. Since we do not expect the
radius of the working surface to be well described by analytical arguments, we
have no reason to expect that the quantity $r_{\rm ws}/\chi^{1/2}$, displayed
in the bottom row of Fig. 4, should be linear.

###  5.3 Evolution of the jet lobes

Having compared the evolution of our jets with simple analytical estimates and
identified regimes in which the momentum flux is not conserved along the jet
axis, we now turn to more quantitative analysis of the energy and mass content
of the jet lobes to better understand the processes driving the morphological
features we have identified.

Throughout this discussion, we will restrict our attention to the ‘jet lobes’
which we define to be all gas cells that have a jet tracer fraction
satisfying:

$f_{\rm J,i}\equiv\frac{m_{\rm J,i}}{m_{\rm i}}>f_{\rm J}^{\rm(thresh)}\;,$
(80)

where we choose the threshold jet fraction to be $f_{\rm
J}^{\rm(thresh)}=10^{-3}$ (Weinberger et al., 2017b; Hardcastle & Krause,
2013). It should be noted, however, that this threshold jet fraction is not an
empirically measurable quantity and, intuitively, we would expect it to be a
complex function of jet power and time, rather than simply a constant. This
work, however, is concerned with better understanding the underlying physical
processes driving the morphology and evolution of jets rather than the precise
determination of the mass and energy content of the jet lobes and we find that
the choice of $f_{\rm J}^{\rm(thresh)}=10^{-3}$ is sufficient for these
purposes.

We begin this analysis by studying the evolution of the total energy and mass
content of the jet lobes, shown in Figs. 6 and 7. As discussed in previous
sections, the jet propagation timescale depends sensitively on the power of
the jet as well as the density of the CGM. In order to make meaningful
comparisons between all jet runs, it is therefore instructive to consider the
evolution of the lobes as a function of the total length of the jet lobes
which acts as a proxy for time, rather than as a function of time itself. To
aid intuition as to the correspondence between jet length and the elapsed
time, on the top $x$-axis of each panel we have plotted blue ticks that are
equally spaced in time. In each plot, the ‘standard’ CGM runs are shown in the
top row and the ‘dense’ CGM runs are shown in the bottom row. Each column
corresponds to the ‘high’, ‘medium’ and ‘low’ power jets, from left to right,
respectively111111Were we to perform the following analysis using data from
both the northern and southern lobes, slight differences in their length
evolution and the location of internal shocks can complicate things. We
therefore restrict the analysis in this section and Section 5.4 to the
northern lobe of each jet and henceforth any reference to the ’jet length’
refers to the length of the northern lobe of the relevant jet and similarly
any quantitative reference to masses and energies are those of the northern
lobe..

Figure 6: The energy composition of the northern lobe of each jet as a
function of the total length of this lobe. The top and bottom rows show the
runs in the ‘standard’ and ‘dense’ CGM cases, respectively, and from left to
right the columns correspond to the ‘high’, ‘medium’ and ‘low’ power jets. In
each panel, the total energy contained within lobe material and the total
energy that has been injected are shown by the thick and thin black lines,
respectively. The orange and purple lines show the kinetic and thermal
components of this total jet lobe energy. Estimates to the thermal energy
coming from CGM and disc material are shown by the dashed and dotted purple
lines, respectively. To aid intuition as to the correspondence between the
relevant jet length and elapsed time, on the top $x$-axis of each panel we
have plotted blue ticks that are equally spaced in time.

####  5.3.1 Lobe energy evolution

In Fig. 6, the total energy contained within the jet lobe and the total energy
that has been injected into the northern lobe of the jet are shown by the
thick and thin black lines, respectively, while the kinetic and thermal
components of this total lobe energy are shown by the orange and purple lines,
respectively121212We wish to emphasise that following analysis should not be
interpreted as predictions for the final partitioning of the injected jet
energy since the volume we define to be the jet lobes can have net
inflow/outflow. Similarly, once turbulence takes hold and disrupts the jet,
our jet tracer definition of ‘lobes’ will not refer to a coherent volume and
so this analysis is better thought of as quantifying the energy content of
jet-dominated material..

We can see that, in most runs, the total energy in the jet lobes remains below
the injected energy throughout their evolution, as would be expected since the
jet inflation does work on the surrounding CGM (and this measure of total
energy does not account for any energy in jet material that does not satisfy
our lobe criterion, although we expect this to have a negligible effect). In
the ‘standard’ CGM runs, when the jet length has reached $3000$ pc, the ratio
of the lobe energy to the injected energy decreases with jet power, ranging
from $0.66$ in the ‘low’ power jet to $0.52$ in the ‘high’ power jet. In the
‘dense’ CGM runs there is not a clear trend, however, as is particularly
evident from the late-time evolution of the ‘low’ power jet, where total
energy in the jet lobes is $1.44$ times the injected energy at the final time
(i.e. when its length has reached $1770\;{\rm pc}$).

The jets in the ‘standard’ CGM case all show similar behaviour: the kinetic
energy component remains dominant throughout (with the average kinetic
fraction ranging from $0.75$ \- $0.85$), the only exception at very early
times when the thermal content of the lobes is high due to initial shock-
heating that occurs as the jet breaks out of the circumnuclear disc. At later
times, even though the thermal energy content never becomes comparable to the
kinetic, as the jets propagate further from the black hole the thermal
component is more significant. The fraction of the lobe energy in thermal
energy is larger for higher power jets and by the time the jet length reaches
$3000\;{\rm pc}$ this amounts to $\sim 0.32$ for the ‘high’ power jet, $\sim
0.26$ for the ‘medium’ power jet and $\sim 0.24$ for the ‘low’ power jet. This
increase in thermal fraction correlates with the strength of the reverse shock
at the head of the jet (see Fig. 5), however other factors such as entrainment
and mixing in the lobes will also be at play here.

Turning now to the runs in the ‘dense’ CGM case, we see that they too are all
dominated by thermal energy during the disc break-out phase and then quickly
transition to kinetic dominance. In contrast to the runs in the ‘standard’
CGM, however, the ‘dense’ CGM runs all rapidly become thermally dominated
again. The time and length-scales over which the transition to thermal
dominance occurs decreases with increasing jet power, ranging from $1.08\;{\rm
Myrs}$ ($Z_{\rm J}=874\;{\rm pc}$) in the ‘low’ jet power case to $0.13\;{\rm
Myrs}$ ($Z_{\rm J}=462\;{\rm pc}$) in the ‘high’ jet power case.

To better understand the behaviour of the thermal energy of the jet lobes we
consider the impact of two distinct processes that primarily influence this
thermal content. Firstly, any mixing or entrainment that occurs, such as in
the backflowing cocoon where shocked jet material mixes with hot CGM, can lead
to an increase in thermal energy. Similarly, instabilities that develop along
the jet beam will also lead to enhanced mixing, as will the entrainment of
disc material. Note that as a result of mixing and entrainment the thermal
energy of the lobes will increase, as will the mass of jet lobes (see Section
5.3.2 below), but only mixed material that satisfies our jet tracer fraction
threshold will contribute to the lobe thermal energy in Fig. 6. The second
channel through which changes to the energy content of the lobes result is
shock heating. This will occur as jet material passes through on-axis shocks
but also in the initial shocks as the jet breaks out of the disc and, where it
exists, in the supersonic turbulence in the lobes.

To better understand how the thermalisation is being driven, for each run we
have provided an estimate (which is likely a lower bound) on the thermal
energy content of the jet lobes that we expect has come from entrained CGM
material, which is indicated by the dashed purple line in Fig. 6. This
quantity has been calculated using our CGM tracer, under the assumption that
the mass associated with this tracer has a specific internal energy equal to
the initial specific internal energy of CGM material (which was initialised
with temperature $T_{\rm CGM}=10^{7}\;{\rm K}$). We have also carried out a
similar analysis using the disc tracer to estimate the contribution of
entrained disc material to the thermal energy of the lobes, under the
assumption that disc material has temperature $T_{\rm disc}=2\times
10^{4}\;{\rm K}$ (again, this is likely a lower bound as we find a non-
negligible amount of disc tracer in the shock-heated cocoon material and some
disc material will have undergone shock heating during the initial relaxation
period). This estimate is shown by the dotted purple line in Fig. 6. These
estimates serve two purposes. Firstly, we can use them to make predictions
about the relative importance of the thermal energy associated with disc and
CGM material in the jet lobes. Secondly, since they are calculated under the
assumption that the material has not undergone shock-heating we can use them
to identify regimes where shock-heating must be important131313It should be
noted, however, that there is degeneracy between these two scenarios (for
example, significant shock-heating of disc material may change its relative
importance). Further, we should point out that it is not possible to perform
similar analysis for the case of ‘pure’ jet material as this would require us
to associate with it a specific internal energy..

We first examine our estimate of thermal energy associated with disc material.
Looking specifically at the ‘high’ power jet in the ‘dense’ CGM case we find
that our disc thermal energy estimate is, on average, $\sim 4$ orders of
magnitude lower than the total thermal energy. In order to explain the total
thermal energy, the disc material would therefore need to have be shock heated
to $\sim 10^{8}\;{\rm K}$. Whilst a non-negligible amount of disc material
ends up in the cocoon in which we do find temperatures that exceed
$10^{8}\;{\rm K}$, the majority of the disc material can be found in the cold
dense sheath surrounding the jet channel (see Fig. 3). For this reason it is
unlikely that the entire thermal energy can be accounted for by entrainment of
disc material, regardless of whether it has undergone shock-heating or not.
Similar arguments can be made for the other runs.

Turning now to our estimate of the contribution to the thermal energy from CGM
material we find that (except at early times) it is much higher in comparison
to our estimate of the disc thermal energy. In both the ‘standard’ and ‘dense’
CGM cases for lower jet powers, the thermal contribution from the CGM material
becomes significant when compared to the total thermal energy of the lobes. At
the final time we find, for these ‘low’ power jets, our CGM thermal energy
estimate is $0.40$ and $0.71$ times the total thermal energy in the ‘standard’
and ‘dense’ cases, respectively. Since our estimates are calculated under the
assumption that the entrained CGM material has not undergone shock-heating we
can infer that, for the thermal energy content to be dominated solely by that
of the entrained CGM material, we must invoke significant heating in the
‘high’ power cases whereas at lower powers, less heating would be needed.

The effects of CGM entrainment are particularly evident in the case of the
‘low’ power jet in the ‘dense’ CGM. When this jet transitions to a state of
thermal energy dominance, the thermal energy of the lobe comes predominantly
from ‘unshocked’ CGM material. This transition occurs as the jet beam starts
to disrupt which is accompanied by the jet length stagnating (as evidenced by
the blue ticks on the $x$-axis). Interestingly, we find that this ultimately
leads to the energy in the lobe of this jet exceeding the total injected
energy.

Figure 7: The mass of the northern lobe of each jet as a function of the total
length of this lobe. The top and bottom rows show the runs in the ‘standard’
and ‘dense’ CGM cases, respectively, and from left to right the columns
correspond to the ‘high’, ‘medium’ and ‘low’ power jets. In each panel, the
total mass of lobe material is shown by the thick black lines. The orange,
purple and grey lines show the the jet, disc and CGM contributions to this
total lobe mass. To aid intuition as to the correspondence between the
relevant jet length and elapsed time, on the top $x$-axis of each panel we
have plotted blue ticks that are equally spaced in time.

####  5.3.2 Lobe mass evolution

We now examine the evolution of the mass of the jet lobes in each of our
simulations, shown by the black line in each of the panels in Fig. 7. Making
use of our passive tracers (see Section 4.2), we have further broken this lobe
mass down into that coming from pure jet, entrained disc and entrained CGM
material which are indicated by the orange, purple and grey lines,
respectively. The general behaviour of the evolution of the total mass of jet
lobes is broadly similar to that of the total energy. Specifically, the mass
of the jet lobes in the ‘standard’ CGM case remains fairly constant, whereas
in the ‘dense’ CGM runs there is an increase, especially in the ‘low’ power
case, as the jet beam is disrupted.

Considering the runs in the ‘standard’ CGM case, it is apparent that the mass
of these jets is entirely dominated by entrained disc material with the time-
averaged fraction of lobe mass in disc material lying in the range $0.85-0.90$
for all of the runs. The mass associated with ‘pure’ jet material is
subdominant at all times (which is unsurprising given the very large
difference in typical densities of these two components). We find that the
mass associated with CGM material is negligible in the initial stages of
evolution, but increases as the jet propagates further from the black hole
when processes such as mixing and entrainment become effective and backflows
develop. At all times, however, the CGM lobe mass content is at least a factor
of $4$ smaller than the contribution from disc material. Comparing the three
different jet powers, it can be seen that the evolution of the mass of each
jet is qualitatively similar, however, with increasing jet power, the time at
which the mass in CGM material exceeds that coming from ‘pure’ jet material
decreases and occurs over shorter length-scales. This is consistent with what
is seen in Fig. 2, specifically that the levels of mixing and turbulence
within the cocoon increase with jet power.

Turning now to the jets in the ‘dense’ CGM case, as mentioned above, the
increase in the total mass of jet lobes is more significant than in the
‘standard’ CGM runs, particularly at later times when the jet propagation
begins to stall in the ‘low’ and ‘medium’ power jets. This gain in mass is
primarily driven by a considerable increase in the mass of CGM material in the
lobes which, while initially subdominant relative to disc material, becomes
the dominant contributor to the lobe mass. This effect is much more pronounced
in the ‘dense’ CGM runs as mixing becomes more effective due to the
development of recollimation shocks and fluid instabilities; in fact the final
percentage contribution to the lobe mass from CGM material lies in the range
$59$-$73$ per cent in the ‘dense’ CGM runs, whereas in the ‘standard’ CGM runs
this same quantity lies in the range $14$-$19$ per cent.

These effects are most pronounced in the case of the ‘low’ power jet where the
stagnation of the jet is associated with a significant increase in the mass of
CGM material in the lobes. This is in line with the previously discussed jet
lobe thermalisation, where contribution from the entrained CGM is more
important in the jets in the ‘dense’ CGM. As is the case for the jets in the
‘standard’ CGM runs, the component of lobe mass associated with ‘pure’ jet
material remains subdominant throughout. We also notice that the timescale
over which the CGM mass fraction becomes dominant with respect to the jet mass
fraction is shorter for higher jet powers ($0.82\;{\rm Myrs}$ in the ‘high’
power jet, $1.68\;{\rm Myrs}$ in the ‘medium’ power jet and $3.84\;{\rm Myrs}$
in the ‘low’ power jet). However, due to the increase in propagation velocity
with jet power, this transition occurs over longer length-scales ($2183\;{\rm
pc}$ in the ‘high’ power jet, $1787\;{\rm pc}$ in the ‘medium’ power jet and
$1567\;{\rm pc}$ in the ‘low’ power jet).

One final point to discuss in Fig. 7 is the fact that the mass of jet material
(material that has, at some point, been in the jet cylinder) in the lobes is
higher in the ‘low’ power jets. We would expect that, to first order, the
total mass of jet material in the simulation scales linearly with time. Since
the plots in Fig. 7 show the mass evolution as a function of jet length it is,
therefore, expected that we would find a higher mass in the slower ‘low’ power
jets. The amount of pure jet material in the simulation will also depend on
the mass of the jet cylinder and, as mentioned in Section 5.2, the cylinder is
initially slightly over-massive until the gas flows in the vicinity of the
black hole are able to be readjust to the jet which takes longer in the ‘low’
power case. This could lead to more massive lobes in the ‘low’ power jets,
although we expect this effect to be relatively small.

###  5.4 Distribution of energy and mass within the jet lobes

Having explored the evolution of the energy and mass content of the jet lobes,
we now look more closely at the distribution of this energy and mass within
the jet lobes at a fixed point in time during evolution of the jet.
Specifically, Figs. 8 and 9 show the normalised energy and mass distributions
of the jet lobe material as a function of distance along the jet axis. In each
plot, the runs in the ‘standard’ CGM are shown in the top row and those in the
‘dense’ CGM are shown in the bottom row. From left to right, each column
corresponds to the ‘high’, ‘medium’ and ‘low’ power jets. The black line in
each panel shows the energy/mass profile. We have also marked the position of
the strongest shock within the jet lobe material with an ‘x’, where the colour
of the marker encodes the Mach number of the shock, using the same colour-
scale as in Fig. 5. As anticipated, we find that in the ‘standard’ CGM runs,
this corresponds to the reverse shock at the head of the jet whereas in the
‘dense’ CGM runs this locates the first recollimation shock. The profile for
each jet was created at the same time as the slices in Figs. 2 and 3 such that
the length of the northern lobe of each jet is $\sim 2\;{\rm kpc}$.

Figure 8: The normalised energy distribution of jet lobe material as a
function of distance along the jet axis. In each plot, the runs in the
‘standard’ CGM are shown in the top row and those in the ‘dense’ CGM are shown
in the bottom row. From left to right, each column corresponds to the ‘high’,
‘medium’ and ‘low’ power jets. The profile for each jet was created at the
same time as the slices in Figs. 2 and 3 meaning that the length of the
northern lobes of the jets are all $\sim 2\;{\rm kpc}$. The black line in each
panel shows the total energy profile. The location of the strongest shock
within the jet lobe material is indicated by an ’x’, where the colour of the
marker encodes the Mach number of the shock, using the same colour-scale as in
Fig. 5. The kinetic and thermal components of this total energy are indicated
by the orange and purple lines, respectively and the purple dashed line
corresponds to our estimate of the thermal energy from the entrained CGM
material. Figure 9: The normalised mass distribution of jet lobe material as a
function of distance along the jet axis. In each plot, the runs in the
‘standard’ CGM are shown in the top row and those in the ‘dense’ CGM are shown
in the bottom row. From left to right, each column corresponds to the ‘high’,
‘medium’ and ‘low’ power jets. The profile for each jet was created at the
same time as the slices in Figs. 2 and 3 meaning that the length of the
northern lobes of the jets are all $\sim 2\;{\rm kpc}$. The black line in each
panel shows the total mass profile. The location of the strongest shock within
the jet lobe material is indicated by an ’x’, where the colour of the marker
encodes the Mach number of the shock, using the same colour-scale as in Fig.
5. The jet, disc and CGM components of this total energy are indicated by the
orange, purple and grey lines, respectively.

####  5.4.1 Energy distribution

In Fig. 8, the total energy distribution is then further broken down into the
kinetic and thermal components, indicated by the orange and purple lines,
respectively. The purple dashed line corresponds to the estimate of the
thermal energy from entrained CGM material (obtained using the same procedure
outlined in Section 5.3.1). We have not included our estimate of the thermal
energy from entrained disc material in Fig. 8 as, in all cases, it was clearly
subdominant along the entire length of the jet.

First, considering the jets in the ‘standard’ CGM runs, we see that the total
jet energy is distributed largely uniformly along the jet out to the reverse
shock. Upstream of this reverse shock, the average kinetic fraction ranges
from $0.87$ in the ‘high’ power jet to $0.92$ in the ‘low’ power which
indicates that these jets are kinetically dominated. Only close to the base of
the jet is the thermal energy content comparable to that of the kinetic energy
which is due to the shock-heating in the disc. Downstream of the reverse
shock, the distribution of total energy in the lobe peaks in the ‘high’ power
jet and we observe an associated peak in the thermal energy distribution,
driven by the thermalisation of material passing through the reverse shock.
Similar behaviour is seen in the ‘medium’ and ‘low’ power jets but the peaks
in the total energy distribution are not as prominent. It should be noted here
that using a jet tracer fraction of $10^{-3}$ to define jet lobe material
means that the majority of the cocoon, including backflowing material, is
categorised as lobe material. The kinetic energy of this backflowing material
will therefore contribute to the profiles in Fig. 8 and explains why in the
‘high’ power jet run, the kinetic energy profile seems to peak downstream of
the strong reverse shock.

The thermal energy of entrained CGM appears largely subdominant for the ‘high’
and ‘medium’ power jet. We find that, upstream of the shock, the average
fraction of thermal energy coming from CGM material is $0.04$ in the ‘high’
power jet and $0.16$ in the ‘medium’ power jet. In the ’low’ power jet,
however, the thermal energy from CGM material does becomes comparable to the
total thermal energy of the jet lobes and we find that, in this case, the
average fraction of thermal energy coming from CGM material in the upstream
gas is $0.43$. Recall however that our estimate of the CGM thermal energy
contribution assumes that this material has not undergone heating. The
observation that, in the ‘high’ power jet, our estimate is significantly lower
than that of the total thermal energy does not mean that the thermal energy of
CGM material is not important here. Rather, it means that if it is to be
important, then it must have undergone shock heating which is indeed more
significant for higher power jets (see Section 5.3).

Turning now to the jets in the ‘dense’ CGM case, we see that as we move along
the jet axis, there is a significant redistribution in the total energy
content, driven by the increasing fraction of energy in thermal form. The
thermal energy then remains dominant out to the head of the jet and we find
that the average fraction of energy in thermal form in the downstream gas
ranges from $0.76$ in the ‘high’ power case to $0.91$ in the ‘low’ power case.
The ‘high’ power jet run shows a significant increase in the contribution of
thermal energy to the total which peaks around the location of the first
recollimation shock and remains approximately constant out to the head of the
jet. Naively, it may be tempting to attribute this increase in thermal energy
distribution to the action of the recollimation shock alone. Whilst material
in the jet channel that passes through this shock will undergo thermalisation,
the slices in Fig. 3 show that there is a significant cocoon of shock-heated
material that extends down almost to the base of the jet. Indeed, the thermal
energy profile begins to increase upstream of the recollimation shock, hence
the on-axis shock heating cannot be the sole actor driving this thermal
profile.

For a significantly higher choice of $f_{\rm J}^{\rm(thresh)}=0.5$ (relative
to our fiducial value of $10^{-3}$) the total lobe energy content of this
‘high’ power jet drops by a factor of $5$ as more of the material in the
turbulent cocoon is no longer characterised as lobe material. Crucially,
however, the value of $(1/E_{\rm lobe})\;{\rm d}E_{\rm lobe,therm}(z)/{\rm
d}z$ at the recollimation shock then drops by a very large factor (of $\sim
18$)141414Note that in the ‘standard’ CGM case the results are much less
sensitive to the adopted choice of $f_{\rm J}^{\rm(thresh)}$ as expected,
given much less developed turbulent cocoon.. This implies that the increase in
thermal energy distribution cannot be attributed solely to the thermalisation
of material passing through this first recollimation shock. Rather, we expect
it to be driven by the presence of hot, off-axis material in the cocoon that
has already undergone shock-heating (potentially by multiple recollimation
shocks)151515One further point to note here is that identifying the peaks in
Mach number in the simulation snapshots may miss any particularly strong
shocks that thermalise on timescales shorter than the interval between
snapshots..

In the ‘high’ and ‘medium’ power jets, the thermal energy distribution
increases significantly close to the base of the jet and then remains fairly
constant from this point out to the head of the jet. This behaviour will be
primarily due to the presence of the hot cocoon that extends down most of the
length of the jets, fed by the series of on-axis recollimation shocks and
associated instabilities which act to thermalise the jet channel material and
divert it off-axis. In the case of the ‘low’ power jet, however, the thermal
energy distribution continues to increase as we move out towards its head.
Comparing this to our estimate of the CGM thermal contribution we see that
very little shock heating is required to explain this total thermal energy
profile, indicating that the trends can be explained by the entrainment of CGM
material that need not have undergone shock-heating.

####  5.4.2 Mass distribution

Fig. 9 shows the distribution of mass within the jet lobe, where the lobe mass
has been split into contributions from disc, CGM and ‘pure’ jet material,
which are indicated by the purple, grey and orange lines, respectively. Note
that for all jets, the mass distribution peaks close to their base. Visual
inspection of Figs. 2 and 3 clearly shows that there is a dense, cold sheath
surrounding the jet channel close to the disc and, indeed, we can see that
this feature in the mass profile can be attributed to entrained disc material.
This cold sheath is slow which explains why there is no analogous peak close
to the base of the jet in the energy profiles shown in Fig. 8. Moving towards
the head of the jet, in all cases the disc component becomes less important,
with disc material primarily residing in the cocoons.

In the ‘standard’ CGM runs, for the ‘high’ and ‘medium’ jet powers, the
contribution from jet material remains subdominant along the entire length of
the jet and as we move towards the head we see that CGM material becomes
important as the turbulence in the cocoon acts to mix CGM material. The ‘low’
power jet remains dominated by disc material for a greater proportion of its
length, with the disc component becoming subdominant relative to the jet
component at $\sim 1130$ pc whereas in the ‘high’ and ‘medium’ power jets the
disc component becomes subdominant relative to the CGM component at $\sim
640\;{\rm pc}$ and $\sim 990\;{\rm pc}$, respectively. This is confirmed by
the slices in Fig. 2 where we can see that the cold sheath surrounding the jet
channel extends furthest along the length of the jet in ‘low’ power case.

In the ‘dense’ CGM case, as expected, the mass coming from CGM material plays
a much more important role here than in the ‘standard’ CGM case, with average
CGM fractions lying in the range $0.70$ to $0.91$ beyond the recollimation
shock. CGM material becomes dominant closer to the base of the jet and remains
so along its entire length, with the contribution from jet material never
being significant at any location along the jet axis. Considering the ‘low’
power jet in particular, there is a broad peak in the mass distribution
towards the head of the jet, which can be associated with mixed CGM material
after disruption of the jet beam as can be clearly seen in the slices of Fig.
3.

##  6 Blandford-Znajek jets: Results

In the previous section we presented simulations of jets with a fixed power
and direction and performed analysis of their launching from a circumnuclear
disc and subsequent propagation through the surrounding CGM. Moving forward
from this, we now examine simulations in which the jet power and direction
evolve according to our full sub-grid model for Blandford-Znajek jets, as
detailed in Sections 2 and 3.

We have carried out a suite of simulations in which we probe the effects of
our choice of model parameters and where we also experiment with different
circumnuclear disc structures to assess the effects of enhanced accretion flow
on the resultant jet morphology and evolution. A full study of the complete
simulation suite is beyond the scope of this work and will be presented in a
companion paper. It is, however, pertinent to employ the insight we have
obtained thus far to make some concrete assessments of how relaxing the
assumptions of fixed jet power and direction will affect the morphology and
evolution of jets. In this section we, therefore, restrict our attention to
three representative simulations which use our full Blandford-Znajek jet
feedback model.

In all of these simulations, the sub-grid $\alpha$-disc is initialised with
mass $M_{\rm d}=10^{3}\;{\rm M_{\odot}}$ and angular momentum direction
$\mathbf{j}_{\rm d}=\mathbf{\hat{z}}$, while the initial black hole mass and
spin are $M_{\rm BH}=10^{6}\;{\rm M_{\odot}}$ and $a=0.2$. This choice of
$\alpha$-disc mass and angular momentum, however, does not guarantee that it
will initially be in equilibrium with respect to torques from the black hole
and inflow from the surroundings. Since the jet properties are sensitive to
the accretion rate onto the black hole and the angular momentum carried by
this material (see equations (2.5) and (39)), we evolve the simulations for
$2\;{\rm Myrs}$ before turning on the jet to allow the sub-grid disc to come