into equilibrium. This reduces the impact of our initial choice of disc mass
and angular momentum direction on the long term properties of the jet. Doing
so introduces further a complication, however, namely that during these
$2\;{\rm Myrs}$ the black hole mass and angular momentum (in particular the
direction of its angular momentum) will also undergo evolution. To mitigate
these problems, we reset the black hole spin magnitude and direction as well
as the accretion disc angular momentum direction to their initial values
immediately before the jet turns on. Note that we do not include these initial
$2\;{\rm Myrs}$ in our analysis and henceforth, any reference to ‘time’ refers
to time elapsed since the jet was turned on.

In the first of the three simulations we present here, the black hole spin is
initially parallel to the $z$-axis and the initial conditions correspond to
our ‘standard’ CGM density. In the second simulation the black hole spin is,
again, initially parallel to the $z$-axis but here we use our ‘dense’ CGM
initial conditions. Finally, in the third simulation, we initialise the black
hole spin direction to be perpendicular to the $z$-axis meaning that the jet
is initially launched into the circumnuclear disc and here we use the initial
conditions with the ‘standard’ density CGM. For clarity, a full accounting of
the parameters of our sub-grid model that differ between runs can be found in
Table 2.

Note that in the derivation of our expression for the Bardeen-Petterson torque
in Section 2.2 we assumed that the angular momenta of the black hole and
$\alpha$-disc were at most ‘slightly misaligned’. In the simulation where the
jet is initially launched into the disc (which we set up intentionally so as
to achieve the maximum effect in terms of torquing) this is evidently no
longer the case. We will perform detailed investigations to determine the
extent to which this assumption affects the torques in our companion paper,
however, we do expect that this assumption does not qualitatively change our
results but may act to increase the spin alignment timescale to some extent.

Table 2: Properties of the simulations that use our full Blandford-Znajek jet model. 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 ‘vertical’ jet launched into the ‘dense’ CGM). The third column corresponds to the initial inclination angle of the jet direction (black hole spin) to the vertical while density of the CGM is listed in the fourth column. Initial jet direction | CGM label | $\big{(}\mathbf{j}_{\rm BH}\cdot\hat{\mathbf{z}}\big{)}_{0}$ | CGM density
---|---|---|---
|  | [∘] | [${\rm g\,cm^{-3}}$]
‘vertical’ | ‘standard’ | $0$ | $10^{-27}$
‘vertical’ | ‘dense’ | $0$ | $10^{-25}$
‘into the disc’ | ‘standard’ | $90$ | $10^{-27}$
Figure 10: Time evolution of various black hole and jet properties for the
three simulations that use our full Blandford-Znajek feedback model. The top-
left panel shows the evolution of the angle between the jet direction and the
$z$-axis. The evolution of the black hole spin parameter, $a$, is shown in the
top-right panel. In the bottom-left panel we show the evolution of the jet
power while the evolution of the mass inflow rate onto the $\alpha$-disc is
shown in the bottom-right panel.

###  6.1 Evolution of the black hole properties

We begin our analysis by describing the evolution of the black hole-accretion
disc system in all three simulations and provide discussion as to which
physical processes are likely to be driving the features that we identify. We
do so with the aid of Fig. 10 which, in the top row shows the evolution of the
angle between the jet direction and the $z$-axis (left-hand panel) and the
evolution of the magnitude of the black hole spin (right-hand panel). In the
bottom row the evolution of the jet power is shown (left-hand panel) as well
as the evolution of the mass inflow rate onto the $\alpha$-disc (right-hand
panel).

In all runs the initial jet power is $\sim 10^{40}\;{\rm erg\,s^{-1}}$, which
is comparable to the ‘low’ power jets161616From equation (39) it is evident
that the jet power depends only on the black hole spin and mass accretion
rate. Since, in all runs, the black hole spin is $0.2$ when the jet turns on,
any difference in the initial jet power comes about due to variations in the
equilibrium mass reached by the $\alpha$-disc during the initial relaxation
period. which were presented in Section 5. Before the launching of the jet
there is sustained mass inflow onto the $\alpha$-disc in all runs171717Note
that evolution before the jet turns on is not shown in Fig. 10.. Upon jet
launching, however, the inflow is very quickly terminated but then resumes at
later times, albeit very intermittently and with varying strengths.
Considering the two ‘vertical’ jet runs, we see that the inflows in the
‘dense’ CGM case resume $\sim 4\;{\rm Myrs}$ later and are less significant
than those found in the ‘standard’ CGM case. We examined radial velocity maps
of the inner regions in these two runs and found that sustained bulk inflows
are less common in the ‘dense’ CGM case due to the higher levels of jet driven
turbulence that prevent the formation of any significant coherent inflow.

Even without external inflow, the black hole continues to accrete from the
$\alpha$-disc. With no mass flowing onto the $\alpha$-disc to replenish the
material accreted by the black hole or entrained by the jet, the mass of the
$\alpha$-disc drops. This, in turn, leads to a decrease in the accretion rate
onto the black hole which, ultimately, causes the jet power to diminish, as is
seen in Fig. 10. In all runs we observe that at later times mass is able to
flow onto the system again, however, this inflow is now very sporadic and
bursty. In the two runs which have a ‘standard’ CGM density, resumption of
inflow occurs after $\sim 1.7\;{\rm Myrs}$ and $\sim 4.0\;{\rm Myrs}$ for the
case where the jet is launched into the disc and the case where it is launched
perpendicular to the disc, respectively. At these times we observe a
significant increase in the power of the relevant jet, as the black hole
accretion rate increases. It is interesting to note that even in our setup
(where the secularly-driven inflow occurs at a moderate rate) the jet power
need not cease entirely in order for the inflow to resume. The jet is
essentially ‘on’ all the time, albeit with its power varying by approximately
one order of magnitude. This illustrates how the jet power is able to
naturally readjust to the properties of the gas in its immediate vicinity and,
in turn, how gas properties are shaped by the jet action (e.g. via shock
heating and entrainment). This suggests that jet duty cycles may not be driven
solely by large-scale inflows (or lack thereof).

Figure 11: Slices in the $x$-$z$ plane (centred on the black hole) for the two
‘vertical’ jets that evolve according to our full Blandford-Znajek feedback
model. The two panels on the left correspond to the run in the ‘standard’ CGM
and the two on the right correspond to the run in the ‘dense’ CGM. For each
run, the first panel is created at an ‘early’ time ($1\;{\rm Myr}$) and the
second at a ‘late’ time ($10\;{\rm Myr}$). The bottom half of each panel shows
the $z$ component of the vorticity field with the colourscale corresponding to
the magnitude of this field and the red and blue colours indicating oppositely
directed vectors. The top half of each panel shows the temperature and onto
this we overlaid the locations of any shocks (where they exist) with the green
colourscale indicating the strength of the shock.

Examining the evolution of spin direction, it is evident that the direction of
the jets whose axes are initialised to be vertical remains approximately
unchanged throughout. This is to be expected as the large-scale gas angular
momentum is well ordered and so it is unlikely that there is significant
accretion of gas that has considerably misaligned angular momenta. In the run
where the jet is initially directed into the disc, however, we see that the
black hole spin (and thus the jet direction) is torqued towards the vertical
as the Bardeen-Petterson effect acts to align it with the total angular
momentum of the black hole–$\alpha$-disc system. The amount of torquing is
significant, with spin direction being changed by $\sim 30^{\circ}$ over the
course of $10\;{\rm Myrs}$. One further point to note here is there is a
distinct knee in the evolution of the spin direction for this misaligned jet
which is coincident with a strong burst of inflow onto the system. This
highlights the fact that, as expected, the rate of alignment of the black hole
angular momentum is also sensitive to the mass flows in the vicinity of the
black hole.

We now turn to look at the evolution of the black hole spin magnitude and,
whilst in all runs it is evident that the spin increases monotonically
throughout, the rate of increase is low. This implies that there is an
interesting interplay between the spin-up torques from accretion of gas from
the corotating $\alpha$-disc and the spin-down torques from the launching of
the Blandford-Znajek jet (we leave detailed characterisation of the relative
importance of these two processes to future work). Furthermore, by
construction, our simulated system can sustain only very moderate inflows onto
the $\alpha$-disc as it is largely in equilibrium (and undergoes only very
modest secular evolution). In more perturbed systems (such as those undergoing
mergers or where large-scale cosmic inflows or cooling flows exist) we would
expect that the black hole accretion rate would be higher which would then
lead to more significant levels of black hole spin-up and higher jet powers.

In this section we have focused on the evolution of the black hole and jet
properties but it is also instructive to understand how the $\alpha$-disc co-
evolves with the black hole. Whilst we save comprehensive analysis for our
future work, it is worth highlighting some key features. As mentioned above,
the $\alpha$-discs are all initialised with a mass of $1000\;{\rm M_{\odot}}$,
but during the first $2\;{\rm Myrs}$ before the jet turns on, in all three
runs the discs settle into quasi-equilibrium with mass of $\sim 400\;{\rm
M_{\odot}}$. Once the jets are launched, the $\alpha$-disc masses decrease
gradually as sporadic inflow events do not supply the discs with sufficient
material to sustain the equilibrium disc masses that were attained when there
was consistent inflow. This leads to a net decrease in the $\alpha$-disc
masses of $\sim 100-200\;{\rm M_{\odot}}$ over the course of the simulations.

When considering the net angular momentum direction of the $\alpha$-disc we
find that the discs in the ‘vertical’ jet cases remain directed along the
$z$-axis, as would be expected since the black hole spins remains
approximately parallel to the $z$-axis and so the disc undergoes no
significant Bardeen-Petterson torque and any inflow from the circumnuclear
disc is largely aligned. In the case of the jet directed into the disc, since
the $\alpha$-disc also initially lies in the plane of the circumnuclear disc,
it is subject to Bardeen-Petterson torques equal and opposite to those felt by
the black hole. These begin to bring its angular momentum into alignment with
the total angular momentum of the black hole-$\alpha$-disc system, albeit
counteracted somewhat by sporadic inflow from the circumnuclear disc. It will
be interesting to probe the effects of the initial $\alpha$-disc configuration
in future work.

###  6.2 Initially vertical jets

We now investigate differences between the jets from Section 5 (in which we
fixed their power and axis) and those whose axes are initially vertical and
evolve according to our full feedback model. For these ‘full’ model runs, we
have highlighted the fact that the spin direction does not change appreciably
(see Fig. 10), hence any observed differences in the jet morphologies with
respect to the ‘fixed’ models are likely due to changes in the jet power
alone. In terms of jet energetics the most relevant runs for comparison are
the ‘low’ power jets in the ‘standard’ and ‘dense’ CGM (i.e. the bottom row of
Figs. 2 and 3).

Figure 12: Slices in the $x$-$z$ plane (centred on the black hole) for the jet
that is initially directed into the circumnuclear disc that evolves according
to our full Blandford-Znajek feedback model. The three smaller panels on the
left show slices of the temperature field, evolving in time with $3\;{\rm
Myr}$ intervals from top to bottom. The main panel shows a slice of the
density field after a further $3\;{\rm Myrs}$ have elapsed. The two inset
plots show slices of the temperature field in the centre and the vorticity
field in the outflow. In the bottom left of each panel the arrow indicates the
direction of the jet (i.e. that of the black hole angular momentum) and this
is labelled with the corresponding inclination angle to the vertical.

In Fig. 11 we show slices in the $x$-$z$ plane (centred on the black hole) for
the two ‘vertical’ jets that evolve according to our ‘full’ model. For each
run we show slices at two different epochs (with the time indicated in the top
right corner). The two columns on the left correspond to the run in the
‘standard’ CGM and the two columns on the right show the run in the ‘dense’
CGM. The bottom half of each panel shows the $z$ component of the vorticity
field with the colourscale corresponding to the magnitude of this field and
the red and blue colours indicate oppositely directed vectors. The top half of
each panel shows the temperature and onto this we have overlaid the locations
of any shocks (where they exist) with the green colourscale indicating the
strength of the shock.

At $1\;{\rm Myr}$ the jet in the ‘standard’ CGM is already undergoing
recollimation and at $10\;{\rm Myrs}$ this process has indeed been completed.
Comparing this to the fixed power runs, it is evident that the recollimation
process occurs over shorter timescales when the jet power is allowed to
decrease in accordance with the accretion rate. The fact that the CGM pressure
in these two runs is the same further reinforces our previous observation
(from comparisons between the ‘fixed’ jets in ‘standard’ and ‘dense’ CGM
cases) that the magnitude of the pressure contrast between the jet lobes and
the surrounding medium is central to the recollimation process. $1\;{\rm Myr}$
after the jet is initially launched, it is still able to drive a bow shock
into the CGM and there is clear evidence for internal shocks in the jet
channel as shells of material that are travelling at different velocities
collide.

At this time, the typical Mach numbers along the jet axis and at the reverse
shock are a factor of $\sim 2.5$ lower than found in the ‘low’ (fixed) power
jet at this time, indicating that the internal shock heating is weaker. This
explains why the channel and lobes of this ‘full’ model jet are colder than
those of the ‘low’ (fixed) power jet. After $10\;{\rm Myrs}$, the ‘full’ model
jet is no longer driving a bow shock into the CGM, however we can still
distinctly see the shock that surrounds the cold jet channel. Interestingly,
there is a plume of material beyond the recollimation shock which has an
enhanced temperature relative to the downstream shocked jet material. This
feature can be attributed to the increase in jet power that occurs due to the
burst of inflow at $\sim 5\;{\rm Myrs}$. This highlights the fact that jet
power variability is likely to enhance levels of mixing, as new jet material
interacts with the previously launched component which will ultimately affect
the morphology and propagation of the jet.

Much of this discussion can also be applied to the ‘full’ jet model in the
‘dense’ CGM but there is one feature specific to this run that is worth
highlighting. Namely, that after $10\;{\rm Myrs}$ the jet beam has been
largely disrupted and mixed with the CGM material. This occurs because the jet
power has dropped by an order of magnitude during the $10\;{\rm Myrs}$ since
its initial launching. Whilst the late-time jet power is too low to impart any
significant heat to the CGM, we see from the vorticity field that the jet is
acting to stir the CGM gas which leads to enhanced levels of turbulence. This
adds weight to the argument that accurate modelling of the full feedback cycle
(including self-regulation) is necessary to improve our understanding of where
and how jet energy is communicated to the surrounding medium and whether jet
heating is localised to the centre of the galaxy or is able to regulate gas
inflows from farther out. We expect that there is no ‘one-size-fits-all’
solution to these problems, but that the dominant processes will be highly
dependent on both the large and small-scale properties of the system in
question. Detailed investigation of these issues will require realistic,
cosmological host galaxy simulations that are evolved with self-regulated jet
feedback.

###  6.3 Jet into the circumnuclear disc

Using the simulation where the jet is initially directed into the
circumnuclear disc, we now briefly consider the effect on the outflow
properties and jet axis evolution when the initial jet direction is not
aligned with the global angular momentum of the system. We perform this
analysis with the aid of Fig. 12 in which the three panels on the left show
slices of the temperature field in the $x$-$z$ plane (again, centred on the
black hole), evolving in time with $3\;{\rm Myr}$ intervals from top to
bottom. The main panel then shows a corresponding slice of the density field
after a further $3\;{\rm Myrs}$ have elapsed with two inset plots showing the
central temperature structure and the vorticity field in the outflow.

From these visualisations, it is immediately obvious that the outflow
structure here is drastically different from all those considered previously
in this work. Here the jet is primarily acting to blow off the surface of the
disc rather than to inflate jet lobes. It is interesting to note, however,
that this outflow still has a somewhat bipolar nature with the majority of the
material being expelled along the $z$-axis, due to the disc pressurisation of
the outflow, as is particularly clear from the temperature slices.

Inspecting the inset vorticity plot in the main panel we see that the amount
of turbulence in the outflow is significantly enhanced beyond that found in
the CGM. We can estimate the properties of the outflow by selecting material
that satisfies our ‘jet lobe condition’ (i.e. $f_{\rm J}>10^{-3}$) and also
has a $z$-velocity greater than $50\;{\rm kms^{-1}}$. After $10\;{\rm Myrs}$
we find that the outflowing material has mass $\sim 7\times 10^{3}\;{\rm
M_{\odot}}$ which corresponds to $~{}0.6$ per cent of the total mass of the
circumnuclear disc. Within this material we also find temperatures and
densities that vary by $\sim 3$ and $\sim 6$ orders of magnitude,
respectively, highlighting the multiphase nature of the outflow. The emergence
of this turbulent multi-phase outflow is interesting for a multitude of
reasons but two in particular stand out. Firstly, this clearly demonstrates
that jets can drive galactic-scale outflows even when directed into a disc,
and secondly, this hints at the possibility of dense clumps forming in the
outflow which could, ultimately, lead to star formation in the outflow itself.
Full investigation of the efficacy of the second scenario is beyond the scope
of this study as such work would require dedicated simulations with realistic
gas cooling (and star formation) in the outflow.

We found in Section 6.1 that in this simulation there is a significant burst
of inflow that feeds the $\alpha$-disc around $\sim 3\;{\rm Myrs}$ after the
initial launching of the jet which results in an increase in the power of the
jet. In the temperature slice created at $4\;{\rm Myrs}$ in Fig. 12 we can see
the result of this increase in jet power in the form of an asymmetric
outflowing shell that has temperatures enhanced above those found in the
extended outflow.

In Section 6.1 we also highlighted the fact that over the course of this
simulation the jet axis is torqued back towards the vertical. To see this
evolution more clearly we have indicated the direction of the jet (i.e. that
of the black hole angular momentum) with an arrow and labelled the inclination
angle in the bottom left of each panel in Fig. 12. At early times the jet is
directed into the disc but due to the effects of sustained accretion from the
$\alpha$-disc the inclination angle decreases, with distinctive asymmetric
features occurring in the outflow. After $10\;{\rm Myrs}$, the jet has begun
to emerge from the disc, as can clearly be seen in the inset temperature plot.

##  7 Discussion

###  7.1 Comparison to previous work

While the majority of AGN jet simulations in the context of galaxy formation
typically consider larger scales and primarily focus on galaxy cluster
environments, it is instructive to compare with these previous works as this
allows us to highlight similarities and differences both in methodology and in
physical outcomes.

Specifically, in this work we aim to improve on existing models of AGN jet
feedback in galaxy-scale simulations by self-consistently predicting the jet
power and direction based on the black hole properties and those of the
surrounding accretion flow. We do so by explicitly following the evolution of
the black hole mass and angular momentum as it accretes from a sub-grid thin
(warped) $\alpha$-disc.

Many jet studies use a simplified prescription for the jet power and direction
and assume both to be fixed (e.g. Bourne & Sijacki, 2020; Bambic & Reynolds,
2019; Ehlert et al., 2018; Bourne & Sijacki, 2017; Weinberger et al., 2017b;
Reynolds et al., 2002). Often this choice is made out of necessity as it is
not always possible to sufficiently resolve the accretion flows. Some works,
however, use more physically motivated models for the jet power and assume
that it is some fraction of the energy associated with the mass accreted by
the black hole. An estimate for the mass accretion rate is then determined
from the gas properties in the simulation and there are a variety of ways of
making such an estimate in the literature; for example, Dubois et al. (2010,
2012) adopt the Bondi-Hoyle-like accretion rate as a measure of the hot gas
accretion rate, while (Yang & Reynolds, 2016b; Li et al., 2017; Li et al.,
2015; Li & Bryan, 2014b) have prescriptions for cold gas accretion via a gas
dropout time. Others use the instantaneous mass accretion rate at the inner
boundary of the simulation domain (Vernaleo & Reynolds, 2006).

More realistic models of black hole spin evolution have begun to be
incorporated into galaxy-scale simulations (Bustamante & Springel, 2019;
Fiacconi et al., 2018; Dubois et al., 2014). More recently, Beckmann et al.
(2019) presented a model in which the jet power was determined by the black
hole spin, making use of the results of GRMHD simulations to obtain a jet
efficiency. There are, however, several differences with respect to our model,
including their accretion prescription which is based on the Bondi-Hoyle flows
and the fact that our model takes into account the back-reaction of the
launching of the jet on the black hole mass.

Vernaleo & Reynolds (2006) performed idealised simulations of jets with fixed
axes that are launched into cluster atmospheres and found that the jet energy
was unable to couple effectively to the ICM. Some works overcome this problem
by including precession (Martizzi et al., 2019; Wang et al., 2019; Ruszkowski
et al., 2017; Li & Bryan, 2014a; Falceta-Gonçalves et al., 2010) and others by
including rapid re-orientation of the jet (Cielo et al., 2018). These
precession and reorientation prescriptions have to be put in ‘by hand’; in our
sub-grid model, however, the jet direction is self-consistently determined by
the black hole spin. Ultimately, this will allow us to assess whether the spin
direction changes we observe in our simulations are sufficient to distribute
the jet energy in such a way that it effectively couples to the ICM. This is
particularly important as recent simulations show that effective communication
of the jet energy to the ICM is possible when simulating more realistic
cluster environments that include turbulence and substructures, without having
to invoke jet precession (e.g. Bourne & Sijacki, 2020, 2017; Morsony et al.,
2010).

In order to accurately evolve the black hole properties and to determine the
mass, energy and momentum injected into the jet, the mass and angular momentum
fluxes in the vicinity of the black hole need to be followed at sufficiently
high resolution. Due to computational constraints this requires adaptively
increasing the resolution close to the black hole. Some works have
successfully implemented techniques to do just this (Angles-Alcazar et al.,
2020; Koudmani et al., 2019; Curtis & Sijacki, 2015) and, further to this,
Bourne & Sijacki (2020, 2017) have successfully demonstrated the use of these
techniques in combination with jet launching.

To accurately follow the inflation of the jet lobes it is essential to ensure
that they are sufficiently resolved. Previous studies have found that to
properly capture the energetics of the jet lobe-ICM interface and the
development of instabilities, additional targeted refinement is often required
(Bourne & Sijacki, 2020; Ehlert et al., 2018; Bourne & Sijacki, 2017;
Weinberger et al., 2017b). Similarly, we also find that maintaining sufficient
resolution in the jet lobes is necessary in order for the dynamics of the jet
to be accurately followed. Additionally, we find that the resolution of the
medium through which the jet propagates also must be sufficiently high too.
Indeed, the evolution of multiphase structures and the resolution required to
accurately capture its dynamics is a seemingly ubiquitous problem in
astrophysics (Bennett & Sijacki, 2020; Fielding et al., 2020; Ji et al., 2019;
Gronke & Oh, 2018; Kim et al., 2017).

It is worth emphasising that there are many methods of injecting jets in
galaxy-scale simulations. In our work, we inject mass, energy and momentum
into a finite mass of material and explicitly follow the inflation of the jet
lobes, as do Bourne & Sijacki (2017); Yang & Reynolds (2016b); Dubois et al.
(2010); Omma et al. (2004). Jet models that are employed in larger scale
cosmological simulations, however, necessarily mimic this lobe inflation phase
and directly inject thermal energy into (off-centre) bubbles (Schaye et al.,
2015; Crain et al., 2015; Sijacki et al., 2007) or impart AGN-driven kinetic
winds (Weinberger et al., 2017a), which bypass the small scale jet
thermalisation physics. Since injecting into a finite mass introduces
numerical mass loading of the jet, some works implement the jet via inflow
boundary conditions (Krause et al., 2012; Morsony et al., 2010; Heinz et al.,
2006).

Finally, the fact that we find such a significant change in outflow morphology
when the jet is directed into the circumnuclear disc highlights the necessity
of accurate modelling of this immediate environment which is typically
neglected in idealised setups and not easily resolvable on cosmological
scales. The quasi-bipolar, turbulent, multi-phase outflow that is launched
from our circumnuclear disc when the jet is directed into it bears some
resemblance to AGN-driven winds (Costa et al., 2020; Faucher-Giguère &
Quataert, 2012; Zubovas & King, 2012; King, 2003) despite the fact that it is,
in fact, jet-driven. This very interesting possibility has not been
investigated much in the theoretical literature.

###  7.2 Caveats

Our sub-grid model is intended for use in galaxy-scale simulations which
(generally) solve the equations of non-relativistic dynamics and in which it
will not be possible to resolve the black hole horizon for the foreseeable
future. This prohibits the entirely self-consistent generation of magnetised
jets and so in our model we chose to specify the flux threading the black hole
horizon using the results of GRMHD simulations from Tchekhovskoy et al.
(2012). As highlighted in Section 2.3, these simulations follow the evolution
of a system in which the accretion disc is in a magnetically arrested state
wherein the magnetic flux on the black hole has reached saturation. Whilst
this is by no means a settled issue, it is reasonable to assume, for the
purposes of parameterising the magnetic flux for this model, that all black
hole accretion discs are, indeed, in a magnetically arrested state (Narayan et
al., 2014).

Our sub-grid disc is modelled as a thin $\alpha$-disc which allows us to
follow the angular momentum flux of accreted material onto the black hole,
meaning that the black hole spin can be accurately evolved. This would not be
possible were we to use an accretion prescription that assumes spherical
symmetry. Knowledge of the black hole spin is required for the Blandford-
Znajek jet model (along with the mass accretion rate), and so our
$\alpha$-disc accretion prescription is what, ultimately, allows us to couple
the accretion flow to the jet properties.

Much theoretical and observational work has gone into the determination of the
properties of black hole accretion flows in systems which are capable of
producing radio jets. A possible picture is that this radio mode feedback
coincides with low Eddington ratios, where electrons and ions decouple into a
two-temperature fluid and mass is accreted through a radiatively inefficient,
thick disc181818Super-Eddington accretion is also likely to lead to the
formation of a thick disc due to the larger optical depths (Abramowicz et al.,
1988). (Narayan & Yi, 1994; Shapiro et al., 1976). In this picture, thin
accretion discs are predominantly associated with quasar mode feedback.

It is not, however, certain whether the existence of a radio-jet necessitates
the presence of a thick accretion disc. Recent optimisations of GRMHD codes
have allowed simulations that resolve the thin disc structure to be performed
(Liska et al., 2020, 2019) and these studies show that thin accretion disc can
indeed launch magnetically dominated jets. Note that some GRMHD simulations
also find that the presence (or lack thereof) of radio jets is highly
dependant on the initial magnetic field configuration which is what primarily
determines whether a large-scale magnetic flux is able to build up on the
black hole (Yuan & Narayan, 2014). Further to this, the quasar-radio dichotomy
primarily stems from observations of X-ray black hole binaries (Fender et al.,
2004) and it is not clear whether this observational picture necessarily
extrapolates to the case of supermassive black holes, with FRII sources being
an obvious counter-example (Merloni & Heinz, 2008; Hardcastle et al., 2006).
Future theoretical effort, in which both thin and thick discs are modelled
self-consistently, will be essential in order to better explore this issue.

In order to understand the workings of our model and in some cases, due to the
additional physical and computational complexity, we neglected several
physical processes in our simulations. Firstly, our simulations were evolved
according to a purely hydrodynamic model which does not include the effects of
magnetic fields. It is expected that magnetic fields act to suppress
instabilities (e.g. Dursi & Pfrommer, 2008; Lyutikov, 2006) and thus may play
a significant role in determining the evolution of the jet-ICM interface. This
has been highlighted in recent jet propagation studies that include magnetic
field evolution (Ehlert et al., 2018; Weinberger et al., 2017b; English et
al., 2016). Neglecting magnetic fields has allowed us to focus specifically on
the workings of our model. It will, however, be important to include such
effects in future work, especially as doing so could provide constraints on
the magnetic flux threading the black hole horizon which could be used to
improve our feedback model.

We also do not model the dynamics of cosmic rays in our simulations. The non-
thermal pressure associated with cosmic ray protons may be significant in the
jet lobes and the excitation and subsequent damping of Alvén waves could play
a significant role in energy transport to the ICM (Brunetti & Jones, 2014;
Enßlin et al., 2007; Brunetti et al., 2004). Recent jet studies have included
the effects of cosmic ray pressure (and transport) (Yang et al., 2019; Ehlert
et al., 2018; Weinberger et al., 2017b; Sijacki et al., 2008). Much work is
also being done to develop improved cosmic ray transport theories (Thomas et
al., 2020; Thomas & Pfrommer, 2019; Jiang & Oh, 2018) and incorporate these
into numerical schemes (Hopkins et al., 2020; Thomas et al., 2020; Dubois et
al., 2019; Pfrommer et al., 2017), meaning that it will be vital to
incorporate cosmic ray physics moving forwards with our model.

Our simulations also do not model any physical viscosity which means that we
are not able to accurately capture the rate at which sound waves are damped
(Berlok et al., 2020; Zweibel et al., 2018; Sijacki & Springel, 2006). It is,
however, important to consider this effect when assessing the relative
importance of the damping of sound waves as a method to communicate the jet
energy to the ICM. Accurate modelling of physical viscosity in weakly
collisional plasma (such as the ICM) is notoriously hard and sensitively
depends on the magnetic field dynamics and topology (Zhuravleva et al., 2019;
Mogavero & Schekochihin, 2014). For similar reasons we have also neglected
thermal conduction (Kannan et al., 2017; Ruszkowski & Oh, 2010) which, in
combination with magnetic fields, is expected to be anisotropic, favoured
along the direction of the magnetic field lines and limited across the
interface that separates the jet from the ICM/CGM due to magnetic draping
(Dursi & Pfrommer, 2008).

##  8 Conclusions

In this work we have developed a self-consistent sub-grid model for black hole
accretion through a (warped) $\alpha$-disc and feedback in the form of a
kinetic Blandford-Znajek jet within the moving mesh code arepo. Our simulation
setup has been chosen to represent an environment analogous to that found at
the centre of radio-loud Seyferts. We considered jets launched by a black hole
located at the centre of a circumnuclear disc. This circumnuclear disc is
itself embedded within a hot CGM and we varied the density of this medium to
consider the effects of increased/decreased pressure on the jet propagation.
We first performed simplified simulations in which we fixed the jet power and
direction, allowing us to more straightforwardly compare our results with
analytic predictions and identify the dominant physical processes at play.
Finally, we examined three exemplary simulation cases in which we employed our
full Blandford-Znajek model, with jets initially perpendicular to or directed
into the plane of the circumnuclear disc.

Our main results regarding the jets with fixed power and direction are:

  * •

Jets that are significantly overpressurised with respect to the surrounding
CGM initially terminate in a strong reverse shock. At later times they tend to
recollimate and this process occurs on a faster time-scale if the pressure
contrast between the jet and the CGM is lower. Recollimation excites
instabilities at the boundary of the jet channel, driving material off-axis
and feeding the surrounding cocoon of diffuse, shock-heated gas. At higher jet
powers, a series of recollimation shocks form along the jet axis as the jet
beam is able to retain significant ram pressure component. At lower powers,
however, the jet beam is disrupted, and in this case the prevailing action of
the jet is to stir the CGM and drive turbulence in the cocoon.

  * •

The majority of jet lobe energy is initially in kinetic form, however, over
the course of the simulations we find that a significant fraction of energy
can thermalise due to shock heating and mixing (primarily with CGM material).
Shock heating is more effective at higher jet powers (and for more
overpressurised jets), whereas at lower powers the thermal energy associated
with mildly shocked (or even unshocked) entrained CGM material is sufficient
to explain the total thermal energy of the jet lobes. Once a large thermal
component builds in the jet lobes this energy can be communicated efficiently
with the surrounding medium.

  * •

The initial evolution of largely overpressurised jets is fairly well described
by the standard analytic models (e.g. Begelman & Cioffi, 1989). Once
recollimation occurs, however, our simulations indicate that a key assumption
inherent in the analytic models is violated: namely, the momentum flux at the
base of the jet is not the same as at the head of the jet and this momentum
flux does evolve with time.

  * •

As jets initially break out from the circumnuclear disc they entrain cold disc
material and in some cases not-insignificant quantities of this disc material
can be present along the entire length of the jet. While energetically this
entrained disc material is largely irrelevant, it dominates the total jet lobe
mass, until (in some cases) sufficient CGM material becomes mixed in, which
occurs largely as a consequence of the excitation of instabilities by the
recollimation process.

Now considering the jets which evolve according to our full Blandford-Znajek
model, the main results are as follows:

  * •

When the power and direction of the AGN jets are determined self-consistently,
the jet morphology and evolution differs significantly from those of jets in
which the power and direction are fixed, even when the jet direction does not
evolve significantly. These differences are substantial enough that the
evolution of these self-regulated jets cannot be predicted using ‘fixed power
and direction’ jet models. Initially, as the jet emerges from the
circumnuclear disc it cuts off gas supply to the black hole and the jet power
consequently drops significantly as the $\alpha$-disc is drained. This leads
to slower jet velocities, weaker shocks and colder material in the jet
channel.

  * •

The jet power is naturally self-regulated: as a consequence of the secular
evolution of the circumnuclear disc, inflowing gas is eventually able to reach
the $\alpha$-disc again causing the jet power to increase significantly.
Hereafter, the inflow is intermittent and bursty, but can sustain a continuous
jet (albeit with varying power) with the supermassive black hole gradually
spinning up.

  * •

For the jet initially launched directly into the circumnuclear disc (as may be
the case if the black hole spin was initially misaligned with respect to the
global angular momentum) the resulting outflow bears little resemblance to the
outcomes previously considered. The jet drives a turbulent, multi-phase,
quasi-bipolar outflow as it blows-off the outer layers of the circumnuclear
disc. Accretion onto the $\alpha$-disc persists and the jet is torqued
efficiently out of the circumnuclear disc. As it emerges, its collimated
outflow propagates into the previously generated multi-phase wind.

Our modelling of Blandford-Znajek jets opens up many new avenues through which
the role of AGN feedback in galaxy formation can be explored. We have shown
here that self-consistent jet power evolution is a natural consequence when
black hole accretion and spin evolution prescriptions are coupled to the jet
feedback. Ultimately, this leads to the energetics and direction of the jet
being determined by, but also acting to shape the properties of the wider
environment. In this work we have explored a scenario in which the central
region of the galaxy is largely in equilibrium, with gas inflows being driven
by modest secular evolution. It will be particularly interesting to apply our
model to more dynamically evolving environments such as major mergers and
cooling flow clusters. This will improve our understanding of the extent to
which self-regulated AGN jet feedback influences the large-scale properties of
these systems and help us determine the dominant processes by which it does
so. We have identified three largely distinct channels through which these
jets shape their surroundings: (i) direct heating via shocks, (ii) mixing and
turbulence, and (iii) multi-phase, galactic-scale outflows in which the
innermost regions of the galactic disc are obliterated and drawn up into the
outflow. A natural next step will be to apply this model in the context of
full cosmological simulations with the ultimate aim of providing stronger
constraints on the relative importance of these interaction channels and
improving our understanding of the role these jets play in shaping the
properties of galaxy populations as a whole.

## Acknowledgements

We would like to thank Martin Haehnelt and Chris Reynolds for helpful comments
on the manuscript. RYT, MAB and DS acknowledge the support from the ERC
starting grant 638707 ‘Black holes and their host galaxies: co-evolution
across cosmic time’ and the STFC. This work was performed using the Cambridge
Service for Data Driven Discovery (CSD3), part of which is operated by the
University of Cambridge Research Computing on behalf of the STFC DiRAC HPC
Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS
capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC
operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure.
This work used the DiRAC@Durham facility managed by the Institute for
Computational Cosmology on behalf of the STFC DiRAC HPC Facility
(www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC
capital grants ST/P002293/1 and ST/R002371/1, Durham University and STFC
operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

## Data availability

The data underlying this article will be shared upon request to the
corresponding author.

## References

  * Abramowicz et al. (1988) Abramowicz M. A., Czerny B., Lasota J. P., Szuszkiewicz E., 1988, ApJ, 332, 646
  * Abramowicz et al. (1995) Abramowicz M. A., Chen X., Kato S., Lasota J.-P., Regev O., 1995, ApJ, 438, L37
  * Angles-Alcazar et al. (2020) Angles-Alcazar D., et al., 2020, arXiv e-prints, p. arXiv:2008.12303
  * Bambic & Reynolds (2019) Bambic C. J., Reynolds C. S., 2019, ApJ, 886, 78
  * Bardeen & Petterson (1975) Bardeen J. M., Petterson J. A., 1975, ApJ, 195, L65
  * Barnes & Hut (1986) Barnes J., Hut P., 1986, Nature, 324, 446
  * Beckmann et al. (2019) Beckmann R. S., et al., 2019, A&A, 631, A60
  * Begelman & Armitage (2014) Begelman M. C., Armitage P. J., 2014, ApJ, 782, L18
  * Begelman & Cioffi (1989) Begelman M. C., Cioffi D. F., 1989, ApJ, 345, L21
  * Bennett & Sijacki (2020) Bennett J. S., Sijacki D., 2020, MNRAS, 499, 597
  * Berlok et al. (2020) Berlok T., Pakmor R., Pfrommer C., 2020, MNRAS, 491, 2919
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
  * Bourne & Sijacki (2017) Bourne M. A., Sijacki D., 2017, MNRAS, 472, 4707
  * Bourne & Sijacki (2020) Bourne M. A., Sijacki D., 2020, arXiv e-prints, p. arXiv:2008.12784
  * Bourne et al. (2019) Bourne M. A., Sijacki D., Puchwein E., 2019, MNRAS, 490, 343
  * Bromberg et al. (2011) Bromberg O., Nakar E., Piran T., Sari R., 2011, ApJ, 740, 100
  * Brunetti & Jones (2014) Brunetti G., Jones T. W., 2014, International Journal of Modern Physics D, 23, 1430007
  * Brunetti et al. (2004) Brunetti G., Blasi P., Cassano R., Gabici S., 2004, MNRAS, 350, 1174
  * Bustamante & Springel (2019) Bustamante S., Springel V., 2019, MNRAS, 490, 4133
  * Cannizzo (2001a) Cannizzo J. K., 2001a, ApJ, 556, 847
  * Cannizzo (2001b) Cannizzo J. K., 2001b, ApJ, 561, L175
  * Cattaneo & Teyssier (2007) Cattaneo A., Teyssier R., 2007, MNRAS, 376, 1547
  * Chatterjee et al. (2019) Chatterjee K., Liska M., Tchekhovskoy A., Markoff S. B., 2019, MNRAS, 490, 2200
  * Chilingarian et al. (2018) Chilingarian I. V., Katkov I. Y., Zolotukhin I. Y., Grishin K. A., Beletsky Y., Boutsia K., Osip D. J., 2018, ApJ, 863, 1
  * Chou et al. (2007) Chou R. C. Y., et al., 2007, The Astrophysical Journal, 670, 116
  * Churazov et al. (2005) Churazov E., Sazonov S., Sunyaev R., Forman W., Jones C., Böhringer H., 2005, MNRAS, 363, L91
  * Cielo et al. (2018) Cielo S., Babul A., Antonuccio-Delogu V., Silk J., Volonteri M., 2018, A&A, 617, A58
  * Clarke et al. (1998) Clarke C. J., Kinney A. L., Pringle J. E., 1998, ApJ, 495, 189
  * Cooke et al. (2016) Cooke K. C., O’Dea C. P., Baum S. A., Tremblay G. R., Cox I. G., Gladders M., 2016, ApJ, 833, 224
  * Costa et al. (2020) Costa T., Pakmor R., Springel V., 2020, MNRAS, 497, 5229
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
  * Croston et al. (2008) Croston J. H., Hardcastle M. J., Kharb P., Kraft R. P., Hota A., 2008, ApJ, 688, 190
  * Curtis & Sijacki (2015) Curtis M., Sijacki D., 2015, MNRAS, 454, 3445
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827
  * Davies et al. (2007) Davies R. I., Müller Sánchez F., Genzel R., Tacconi L. J., Hicks E. K. S., Friedrich S., Sternberg A., 2007, ApJ, 671, 1388
  * Donahue et al. (2010) Donahue M., et al., 2010, ApJ, 715, 881
  * Dotti et al. (2013) Dotti M., Colpi M., Pallini S., Perego A., Volonteri M., 2013, ApJ, 762, 68
  * Downes & Solomon (1998) Downes D., Solomon P. M., 1998, ApJ, 507, 615
  * Dubois et al. (2010) Dubois Y., Devriendt J., Slyz A., Teyssier R., 2010, MNRAS, 409, 985
  * Dubois et al. (2012) Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012, MNRAS, 420, 2662
  * Dubois et al. (2014) Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., 2014, MNRAS, 440, 2333
  * Dubois et al. (2019) Dubois Y., Commerçon B., Marcowith A. r., Brahimi L., 2019, A&A, 631, A121
  * Dursi & Pfrommer (2008) Dursi L. J., Pfrommer C., 2008, ApJ, 677, 993
  * Ehlert et al. (2018) Ehlert K., Weinberger R., Pfrommer C., Pakmor R., Springel V., 2018, MNRAS, 481, 2878
  * English et al. (2016) English W., Hardcastle M. J., Krause M. G. H., 2016, MNRAS, 461, 2025
  * Enßlin et al. (2007) Enßlin T. A., Pfrommer C., Springel V., Jubelgas M., 2007, A&A, 473, 41
  * Fabian (1994) Fabian A. C., 1994, ARA&A, 32, 277
  * Falceta-Gonçalves et al. (2010) Falceta-Gonçalves D., Caproni A., Abraham Z., Teixeira D. M., de Gouveia Dal Pino E. M., 2010, ApJ, 713, L74
  * Falle (1991) Falle S. A. E. G., 1991, MNRAS, 250, 581
  * Faucher-Giguère & Quataert (2012) Faucher-Giguère C.-A., Quataert E., 2012, MNRAS, 425, 605
  * Fender et al. (2004) Fender R. P., Belloni T. M., Gallo E., 2004, MNRAS, 355, 1105
  * Fiacconi et al. (2018) Fiacconi D., Sijacki D., Pringle J. E., 2018, MNRAS, 477, 3807
  * Fielding et al. (2020) Fielding D. B., Ostriker E. C., Bryan G. L., Jermyn A. S., 2020, ApJ, 894, L24
  * Fogarty et al. (2019) Fogarty K., et al., 2019, ApJ, 879, 103
  * Gallimore et al. (2006) Gallimore J. F., Axon D. J., O’Dea C. P., Baum S. A., Pedlar A., 2006, AJ, 132, 546
  * García-Burillo & Combes (2012) García-Burillo S., Combes F., 2012, in Journal of Physics Conference Series. p. 012050 (arXiv:1205.0758), doi:10.1088/1742-6596/372/1/012050
  * Gaspari et al. (2011) Gaspari M., Melioli C., Brighenti F., D’Ercole A., 2011, MNRAS, 411, 349
  * Goodman (2003) Goodman J., 2003, MNRAS, 339, 937
  * Greene et al. (2019) Greene J. E., Strader J., Ho L. C., 2019, arXiv e-prints, p. arXiv:1911.09678
  * Gronke & Oh (2018) Gronke M., Oh S. P., 2018, MNRAS, 480, L111
  * Hardcastle & Krause (2013) Hardcastle M. J., Krause M. G. H., 2013, MNRAS, 430, 174
  * Hardcastle et al. (2006) Hardcastle M. J., Evans D. A., Croston J. H., 2006, MNRAS, 370, 1893
  * Harrison et al. (2018) Harrison R., Gottlieb O., Nakar E., 2018, MNRAS, 477, 2128
  * Hawley & Krolik (2006) Hawley J. F., Krolik J. H., 2006, ApJ, 641, 103
  * Heinz et al. (2006) Heinz S., Brüggen M., Young A., Levesque E., 2006, MNRAS, 373, L65
  * Henden et al. (2018) Henden N. A., Puchwein E., Shen S., Sijacki D., 2018, MNRAS, 479, 5385
  * Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  * Hicks et al. (2013) Hicks E. K. S., Davies R. I., Maciejewski W., Emsellem E., Malkan M. A., Dumas G., Müller-Sánchez F., Rivers A., 2013, ApJ, 768, 107
  * Hopkins et al. (2020) Hopkins P. F., et al., 2020, MNRAS, 492, 3465
  * Hota & Saikia (2006) Hota A., Saikia D. J., 2006, MNRAS, 371, 945
  * Imanishi & Nakanishi (2013) Imanishi M., Nakanishi K., 2013, AJ, 146, 47
  * Ji et al. (2019) Ji S., Oh S. P., Masterson P., 2019, MNRAS, 487, 737
  * Jiang & Oh (2018) Jiang Y.-F., Oh S. P., 2018, ApJ, 854, 5
  * Kaiser & Alexander (1997) Kaiser C. R., Alexander P., 1997, MNRAS, 286, 215
  * Kannan et al. (2017) Kannan R., Vogelsberger M., Pfrommer C., Weinberger R., Springel V., Hernquist L., Puchwein E., Pakmor R., 2017, ApJ, 837, L18
  * Kim et al. (2017) Kim C.-G., Ostriker E. C., Raileanu R., 2017, ApJ, 834, 25
  * King (2003) King A., 2003, ApJ, 596, L27
  * King et al. (2005) King A. R., Lubow S. H., Ogilvie G. I., Pringle J. E., 2005, MNRAS, 363, 49
  * King et al. (2007) King A. R., Pringle J. E., Livio M., 2007, MNRAS, 376, 1740
  * Kinney et al. (2000) Kinney A. L., Schmitt H. R., Clarke C. J., Pringle J. E., Ulvestad J. S., Antonucci R. R. J., 2000, ApJ, 537, 152
  * Komissarov (2001) Komissarov S. S., 2001, MNRAS, 326, L41
  * Koudmani et al. (2019) Koudmani S., Sijacki D., Bourne M. A., Smith M. C., 2019, MNRAS, 484, 2047
  * Krause et al. (2012) Krause M., Alexander P., Riley J., Hopton D., 2012, MNRAS, 427, 3196
  * Ledlow et al. (2001) Ledlow M. J., Owen F. N., Yun M. S., Hill J. M., 2001, ApJ, 552, 120
  * Lense & Thirring (1918) Lense J., Thirring H., 1918, Physikalische Zeitschrift, 19, 156
  * Li & Bryan (2014a) Li Y., Bryan G. L., 2014a, ApJ, 789, 54
  * Li & Bryan (2014b) Li Y., Bryan G. L., 2014b, ApJ, 789, 153
  * Li et al. (2015) Li Y., Bryan G. L., Ruszkowski M., Voit G. M., O’Shea B. W., Donahue M., 2015, ApJ, 811, 73
  * Li et al. (2017) Li Y., Ruszkowski M., Bryan G. L., 2017, ApJ, 847, 106
  * Liska et al. (2018) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S., 2018, MNRAS, 474, L81
  * Liska et al. (2019) Liska M., Tchekhovskoy A., Ingram A., van der Klis M., 2019, MNRAS, 487, 550
  * Liska et al. (2020) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S. B., Van Moer M., 2020, MNRAS,
  * Lodato & Pringle (2006) Lodato G., Pringle J. E., 2006, MNRAS, 368, 1196
  * Lyutikov (2006) Lyutikov M., 2006, MNRAS, 373, 73
  * Mao et al. (2015) Mao M. Y., et al., 2015, MNRAS, 446, 4176
  * Martin et al. (2007) Martin R. G., Pringle J. E., Tout C. A., 2007, MNRAS, 381, 1617
  * Martizzi et al. (2019) Martizzi D., Quataert E., Faucher-Giguère C.-A., Fielding D., 2019, MNRAS, 483, 2465
  * McCarthy et al. (2017) McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2017, MNRAS, 465, 2936
  * McKinney (2006) McKinney J. C., 2006, MNRAS, 368, 1561
  * McKinney & Gammie (2004) McKinney J. C., Gammie C. F., 2004, ApJ, 611, 977
  * McKinney et al. (2012) McKinney J. C., Tchekhovskoy A., Bland ford R. D., 2012, MNRAS, 423, 3083
  * McNamara & Nulsen (2007) McNamara B. R., Nulsen P. E. J., 2007, ARA&A, 45, 117
  * McNamara & Nulsen (2012) McNamara B. R., Nulsen P. E. J., 2012, New Journal of Physics, 14, 055023
  * Medling et al. (2014) Medling A. M., et al., 2014, ApJ, 784, 70
  * Merloni & Heinz (2008) Merloni A., Heinz S., 2008, MNRAS, 388, 1011
  * Mezcua et al. (2018) Mezcua M., Civano F., Marchesi S., Suh H., Fabbiano G., Volonteri M., 2018, MNRAS, 478, 2576
  * Mezcua et al. (2019) Mezcua M., Suh H., Civano F., 2019, MNRAS, 488, 685
  * Michel (1973) Michel F. C., 1973, ApJ, 180, L133
  * Mingo et al. (2011) Mingo B., Hardcastle M. J., Croston J. H., Evans D. A., Hota A., Kharb P., Kraft R. P., 2011, ApJ, 731, 21
  * Mingo et al. (2012) Mingo B., Hardcastle M. J., Croston J. H., Evans D. A., Kharb P., Kraft R. P., Lenc E., 2012, ApJ, 758, 95
  * Mogavero & Schekochihin (2014) Mogavero F., Schekochihin A. A., 2014, MNRAS, 440, 3226
  * Morganti et al. (1999) Morganti R., Killeen N. E. B., Ekers R. D., Oosterloo T. A., 1999, MNRAS, 307, 750
  * Morsony et al. (2010) Morsony B. J., Heinz S., Brüggen M., Ruszkowski M., 2010, MNRAS, 407, 1277
  * Nagar & Wilson (1999) Nagar N. M., Wilson A. S., 1999, ApJ, 516, 97
  * Narayan & Medvedev (2001) Narayan R., Medvedev M. V., 2001, ApJ, 562, L129
  * Narayan & Yi (1994) Narayan R., Yi I., 1994, ApJ, 428, L13
  * Narayan & Yi (1995) Narayan R., Yi I., 1995, ApJ, 444, 231
  * Narayan et al. (2014) Narayan R., McClintock J. E., Tchekhovskoy A., 2014, Energy Extraction from Spinning Black Holes Via Relativistic Jets, doi:10.1007/978-3-319-06349-2_25.
  * Novikov & Thorne (1973) Novikov I. D., Thorne K. S., 1973, in Black Holes (Les Astres Occlus). pp 343–450
  * Omma et al. (2004) Omma H., Binney J., Bryan G., Slyz A., 2004, MNRAS, 348, 1105
  * Ostriker et al. (2010) Ostriker J. P., Choi E., Ciotti L., Novak G. S., Proga D., 2010, ApJ, 722, 642
  * Pakmor et al. (2016) Pakmor R., Springel V., Bauer A., Mocz P., Munoz D. J., Ohlmann S. T., Schaal K., Zhu C., 2016, MNRAS, 455, 1134
  * Penna et al. (2010) Penna R. F., McKinney J. C., Narayan R., Tchekhovskoy A., Shafee R., McClintock J. E., 2010, MNRAS, 408, 752
  * Penna et al. (2013) Penna R. F., Narayan R., Sądowski A. e., 2013, MNRAS, 436, 3741
  * Penny et al. (2018) Penny S. J., et al., 2018, MNRAS, 476, 979
  * Perego et al. (2009) Perego A., Dotti M., Colpi M., Volonteri M., 2009, MNRAS, 399, 2249
  * Perucho & Martí (2007) Perucho M., Martí J. M., 2007, MNRAS, 382, 526
  * Peterson & Fabian (2006) Peterson J. R., Fabian A. C., 2006, Phys. Rep., 427, 1
  * Pfrommer et al. (2017) Pfrommer C., Pakmor R., Schaal K., Simpson C. M., Springel V., 2017, MNRAS, 465, 4500
  * Pringle (1992) Pringle J. E., 1992, MNRAS, 258, 811
  * Pringle (1996) Pringle J. E., 1996, MNRAS, 281, 357
  * Pringle (1997) Pringle J. E., 1997, MNRAS, 292, 136
  * Reynolds et al. (2002) Reynolds C. S., Heinz S., Begelman M. C., 2002, MNRAS, 332, 271
  * Reynolds et al. (2015) Reynolds C. S., Balbus S. A., Schekochihin A. A., 2015, ApJ, 815, 41
  * Russell et al. (2017) Russell H. R., et al., 2017, ApJ, 836, 130
  * Russell et al. (2019) Russell H. R., et al., 2019, MNRAS, 490, 3025
  * Ruszkowski & Oh (2010) Ruszkowski M., Oh S. P., 2010, ApJ, 713, 1332
  * Ruszkowski et al. (2017) Ruszkowski M., Yang H. Y. K., Reynolds C. S., 2017, ApJ, 844, 13
  * Schaal & Springel (2015) Schaal K., Springel V., 2015, MNRAS, 446, 3992
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Schinnerer et al. (1999) Schinnerer E., Eckart A., Tacconi L. J., 1999, The Astrophysical Journal, 524, L5
  * Schmitt et al. (1997) Schmitt H. R., Kinney A. L., Storchi-Bergmann T., Antonucci Robert 1997, ApJ, 477, 623
  * Schreiber et al. (2003) Schreiber M. R., Hameury J. M., Lasota J. P., 2003, A&A, 410, 239
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Shapiro et al. (1976) Shapiro S. L., Lightman A. P., Eardley D. M., 1976, ApJ, 204, 187
  * Sijacki & Springel (2006) Sijacki D., Springel V., 2006, MNRAS, 371, 1025
  * Sijacki et al. (2007) Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877
  * Sijacki et al. (2008) Sijacki D., Pfrommer C., Springel V., Enßlin T. A., 2008, MNRAS, 387, 1403
  * Sijacki et al. (2015) Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575
  * Silk (2017) Silk J., 2017, ApJ, 839, L13
  * Sądowski & Narayan (2015) Sądowski A., Narayan R., 2015, MNRAS, 453, 3213
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
  * Tanabe & Nagataki (2008) Tanabe K., Nagataki S., 2008, Phys. Rev. D, 78, 024004
  * Tchekhovskoy et al. (2010) Tchekhovskoy A., Narayan R., McKinney J. C., 2010, ApJ, 711, 50
  * Tchekhovskoy et al. (2011) Tchekhovskoy A., Narayan R., McKinney J. C., 2011, MNRAS, 418, L79
  * Tchekhovskoy et al. (2012) Tchekhovskoy A., McKinney J. C., Narayan R., 2012, in Journal of Physics Conference Series. p. 012040 (arXiv:1202.2864), doi:10.1088/1742-6596/372/1/012040
  * Thomas & Pfrommer (2019) Thomas T., Pfrommer C., 2019, MNRAS, 485, 2977
  * Thomas et al. (2020) Thomas T., Pfrommer C., Pakmor R., 2020, arXiv e-prints, p. arXiv:2010.11960
  * Thorne (1974) Thorne K. S., 1974, ApJ, 191, 507
  * Vernaleo & Reynolds (2006) Vernaleo J. C., Reynolds C. S., 2006, ApJ, 645, 83
  * Voigt & Fabian (2004) Voigt L. M., Fabian A. C., 2004, MNRAS, 347, 1130
  * Wang et al. (2019) Wang C., Li Y., Ruszkowski M., 2019, MNRAS, 482, 3576
  * Weinberger et al. (2017a) Weinberger R., et al., 2017a, MNRAS, 465, 3291
  * Weinberger et al. (2017b) Weinberger R., Ehlert K., Pfrommer C., Pakmor R., Springel V., 2017b, MNRAS, 470, 4530
  * Weinberger et al. (2018) Weinberger R., et al., 2018, MNRAS, 479, 4056
  * Weinberger et al. (2020) Weinberger R., Springel V., Pakmor R., 2020, ApJS, 248, 32
  * Williams et al. (2017) Williams D. R. A., et al., 2017, MNRAS, 472, 3842
  * Yang & Reynolds (2016a) Yang H. Y. K., Reynolds C. S., 2016a, ApJ, 818, 181
  * Yang & Reynolds (2016b) Yang H. Y. K., Reynolds C. S., 2016b, ApJ, 829, 90
  * Yang et al. (2019) Yang H. Y. K., Gaspari M., Marlow C., 2019, ApJ, 871, 6
  * Yuan & Narayan (2014) Yuan F., Narayan R., 2014, ARA&A, 52, 529
  * Zakamska & Narayan (2003) Zakamska N. L., Narayan R., 2003, ApJ, 582, 162
  * Zhuravleva et al. (2019) Zhuravleva I., Churazov E., Schekochihin A. A., Allen S. W., Vikhlinin A., Werner N., 2019, Nature Astronomy, 3, 832
  * Zubovas & King (2012) Zubovas K., King A., 2012, ApJ, 745, L34
  * Zweibel et al. (2018) Zweibel E. G., Mirnov V. V., Ruszkowski M., Reynolds C. S., Yang H. Y. K., Fabian A. C., 2018, ApJ, 858, 5

##  Appendix A Spin dependent parameters

The spin dependent ISCO radius is given by

$R_{\rm ISCO}(a)=\Lambda(a)\frac{G\,M_{\rm BH}}{c^{2}}\,,$ (81)

where

$\Lambda(a)=3+Z_{2}(a)\mp\sqrt{(3-Z_{1}(a))(3+Z_{1}(a)+2Z_{2}(a))}\,.$ (82)

In this expression the upper and lower subtraction and addition operators
differentiate between prograde and retrograde motion cases, correspondingly.
$Z_{1}(a)$ and $Z_{2}(a)$ are given by

$Z_{1}(a)=1+(1-a^{2})^{1/3}\big{[}(1+a)^{1/3}+(1-a)^{1/3}\big{]}\,,$ (83)

and

$Z_{2}(a)=\sqrt{3a^{2}+Z_{1}^{2}(a)}.$ (84)

The spin dependent radiative efficiency is then

$\epsilon_{\rm r}(a)=1-\sqrt{1-\frac{2}{3\Lambda(a)}}\,.$ (85)

Finally, the specific angular momentum at the ISCO is

$L_{\rm ISCO}(a)=\pm\frac{GM_{\rm BH}}{c\Lambda(a)}\frac{\Lambda(a)^{2}\mp
2a\sqrt{\Lambda(a)}+a^{2}}{\sqrt{\Lambda(a)-3\pm 2a/\sqrt{\Lambda(a)}}}\,,$
(86)

where, again, the upper and lower addition and subtraction operators
differentiate between prograde and retrograde motion.

##  Appendix B The Blandford-Znajek energy and angular momentum fluxes

Throughout this section we use comma and semicolon notation for partial and
covariant derivatives. For clarity, we use units such that $G=M_{\rm BH}=c=1$
in the initial calculations but reintroduce physical units at the end.
Throughout the derivation, we work in Kerr-Schild (horizon-penetrating)
coordinates. Further, we assume that the fields are stationary, axisymmetric,
symmetric about $\theta=\pi/2$ and that the black hole magnetosphere is force-
free.

The quantities we wish to determine are the radial fluxes of energy and
angular momentum across the black hole horizon, as measured by a stationary
observer at infinity

$\displaystyle\dot{E}_{\rm BZ}=-\int_{r=r_{\rm
H}}T^{r}{}_{t}\sqrt{-g}\;\mathrm{d}\theta\mathrm{d}\phi\,$ (87)
$\displaystyle\dot{J}_{\rm BZ}=-\int_{r=r_{\rm
H}}T^{r}{}_{\phi}\sqrt{-g}\;\mathrm{d}\theta\mathrm{d}\phi\,,$ (88)

where $T^{\mu\nu}$ is the stress-energy tensor, $r_{\rm H}=r_{\rm
g}(1+\sqrt{1-a^{2}})$ is the radial coordinate of the black hole horizon and
$g$ is the determinant of $g^{\mu\nu}$, the metric tensor.

Under the force-free condition $T^{\mu\nu}$ reduces to the electromagnetic
stress-energy tensor, $T_{\rm EM}^{\mu\nu}$, which is given by

$T_{\rm
EM}^{\mu\nu}=F^{\mu\gamma}F^{\nu}{}_{\gamma}-\frac{1}{4}g^{\mu\nu}F^{\alpha\beta}F_{\alpha\beta}\;,$
(89)

where $F_{\mu\nu}$ is the Faraday tensor191919We have absorbed a factor of
$\frac{1}{\sqrt{4\pi}}$ into the definition of the Faraday tensor. which, in
terms of the electromagnetic four-potential, $A^{\mu}$, is

$F_{\mu\nu}=A_{\nu,\mu}-A_{\mu,\nu}\;,$ (90)

Henceforth we will drop the subscript ‘EM’ in the electromagnetic stress-
energy tensor.

The assumption of ideal MHD gives the constraint

${}^{*}F^{\mu\nu}F_{\mu\nu}=0\,,$ (91)

where ${}^{*}F^{\mu\nu}$ is the dual of the Faraday tensor. It then follows
that

$A_{\phi,\theta}A_{t,r}=A_{t,\theta}A_{\phi,r}\,,$ (92)

and $A_{t}$ is a function of $A_{\phi}$ and therefore constant on magnetic
surfaces (surfaces defined by $A_{\phi}=\text{const}$). Thus, we define

$-\omega(r,\theta)\equiv\frac{A_{t,\theta}}{A_{\phi,\theta}}=\frac{A_{t,r}}{A_{\phi,r}}\,,$
(93)

such that $\omega(A_{\phi})$ is a function that is also constant on magnetic
surfaces.

Imposing stationarity and axisymmetry then gives $F_{\phi t}=F_{t\phi}=0$. The
remaining non-zero components of the Faraday tensor can then be written in
terms of the three variables $A_{\phi}$, $\omega$ and $B^{\phi}$ as

$\displaystyle F_{r\phi}=-F_{\phi r}=A_{\phi,r}\,,$ $\displaystyle
F_{\theta\phi}=-F_{\phi\theta}=A_{\phi,\theta}\,,$ $\displaystyle
F_{tr}=-F_{rt}=\omega A_{\phi,r}\,,$ $\displaystyle F_{t\theta}=-F_{\theta
t}=\omega A_{\phi,\theta}\,,$ $\displaystyle F_{r\theta}=-F_{\theta
r}=\sqrt{-g}B^{\phi}\,,$ (94)

and we also identify

$\displaystyle B^{r}=\frac{1}{\sqrt{-g}}A_{\phi,\theta}\,,$ $\displaystyle
B^{\theta}=-\frac{1}{\sqrt{-g}}A_{\phi,r}\,.$ (95)

The relevant components of the stress-energy tensor, evaluated on the black
hole horizon, are then

$\displaystyle
T^{r}{}_{t}=-2\bigg{(}\frac{A_{\phi,\theta}}{\sqrt{-g}}\bigg{)}^{2}r_{\rm
H}\omega(\omega-\Omega_{\rm H})\sin^{2}\theta\,,$ (96) $\displaystyle
T^{r}{}_{\phi}=\frac{1}{\omega}\;T^{r}{}_{t}\;.$ (97)

The evaluation of equations (87) and (88) therefore require the determination
of the fields $A_{\phi}$, $\omega$ and $B^{\phi}$. This is done by solving the
equation of motion

$T^{\mu\nu}_{\rm EM\,;\nu}=0\;.$ (98)

The resulting equations are highly non-linear and we hence proceed, as in
Blandford & Znajek (1977), by perturbing the known split-monopole solution for
the non-spinning case (Schwarzschild solution) and take quantities to some
order of in $a$.

For this $a=0$ case, the equation of motion reduces to

$\mathcal{L}A_{\phi}^{(0)}=0\,,$ (99)

where we use the superscript $(0)$ to indicate solutions to this non-spinning
case and $\mathcal{L}$ is the linear operator

$\mathcal{L}\equiv\frac{1}{\sin\theta}\frac{\partial}{\partial{r}}\bigg{[}\bigg{(}1-\frac{2}{r}\bigg{)}\frac{\partial}{\partial{r}}\bigg{]}+\frac{1}{r^{2}}\frac{\partial}{\partial{\theta}}\bigg{[}\frac{1}{\sin\theta}\frac{\partial}{\partial{\theta}}\bigg{]}\,.$
(100)

The solutions for a split-monopole field geometry are

$\displaystyle A_{\phi}^{(0)}(r,\theta)=-C\cos\theta\,,$ (101)
$\displaystyle\omega^{(0)}(r,\theta)=0\,,$ (102) $\displaystyle
B^{\phi(0)}(r,\theta)=0\,,$ (103)

where $C$ is a constant.

Perturbing this to the next highest order in $a$, the fields take the form

$\displaystyle
A_{\phi}(r,\theta)=-C\cos\theta+A_{\phi}^{(2)}(r,\theta)a^{2}+\mathcal{O}(a^{4})\,,$
(104)
$\displaystyle\omega(r,\theta)=\omega^{(1)}(r,\theta)a+\mathcal{O}(a^{3})\,,$
(105) $\displaystyle
B^{\phi}(r,\theta)=B^{\phi(1)}(r,\theta)a+\mathcal{O}(a^{3})\,.$ (106)

$A_{\phi}$ is even in $a$ since the radial magnetic field geometry should not
depend on the direction of the rotation. Symmetry arguments further constrain
$\omega$ and $B^{\phi}$ to be odd functions of $a$.
To $\mathcal{O}(a^{2})$ the equation of motion is

$\mathcal{L}A_{\phi}^{(2)}=S(r,\theta)\,,$ (107)

where $S(r,\theta)$ is a source function that depends on $\omega^{(1)}$ and
$B^{\phi(1)}$. Blandford & Znajek (1977) proceed by requiring that the fields
are finite on the horizon and then match this solution to the flat-space
solution at infinity (Michel, 1973).

Imposing these conditions specifies $\omega^{(1)}$ and $B^{\phi(1)}$, and
$A_{\phi}^{(2)}$ can thus be found by inverting equation (107) using the known
Green’s function for the linear operator. Altogether this gives

$\displaystyle
A_{\phi}(r,\theta)=-C\cos\theta+Cf(r)\cos\theta\sin^{2}\theta\,a^{2}+\mathcal{O}(a^{4})\,,$
(108) $\displaystyle\omega(r,\theta)=\frac{1}{8}a+\mathcal{O}(a^{3})\,,$ (109)
$\displaystyle
B^{\phi}(r,\theta)=-\frac{C}{8r^{2}}\bigg{(}1+\frac{4}{r}\bigg{)}a+\mathcal{O}(a^{3})\,,$
(110)

where

$\displaystyle f(r)=$
$\displaystyle\Bigg{[}Li_{2}\bigg{(}\frac{2}{r}\bigg{)}-\ln\bigg{(}1-\frac{2}{r}\bigg{)}\ln\bigg{(}\frac{r}{2}\bigg{)}\Bigg{]}\frac{r^{2}(2r-3)}{8}$
(111)
$\displaystyle+\frac{1+3r-6r^{2}}{12}\ln\bigg{(}\frac{r}{2}\bigg{)}+\frac{11}{72}+\frac{1}{3r}+\frac{r}{2}-\frac{r^{2}}{2}\,,$
(112)

and $Li_{2}(x)$ is the second polylogarithm function

$Li_{2}(x)=-\int_{0}^{1}\frac{\ln(1-tx)}{t}\mathrm{d}t\,.$ (113)

Since the flux integrals only require knowledge of the horizon fields, it is
important to note that the value of $f$ evaluated at $r=r_{\rm H}$

$f(r_{\rm H})=\frac{6\pi^{2}-49}{72}\,.$ (114)

The first correction to the outward fluxes of electromagnetic energy and
angular momentum are then found by inserting these results into equations (87)
and (88), and evaluating the resulting integrals to the lowest order in $a$

$\displaystyle\dot{E}_{\rm BZ}$
$\displaystyle=8\pi\int_{0}^{\frac{\pi}{2}}\frac{1}{\sqrt{-g}}A_{\phi,\theta}^{2}\omega(\Omega_{\rm
H}-\omega)r_{\rm H}\sin^{2}\theta\,\mathrm{d}\theta$
$\displaystyle=\frac{\pi}{24}a^{2}C^{2}+\mathcal{O}(a^{4})\,,$ (115)

and

$\displaystyle\dot{J}_{\rm BZ}$
$\displaystyle=8\pi\int_{0}^{\frac{\pi}{2}}\frac{1}{\sqrt{-g}}A_{\phi,\theta}^{2}(\Omega_{H}-\omega)\sin^{2}\theta\,r_{\rm
H}\mathrm{d}\theta$ $\displaystyle=\frac{\pi}{3}aC^{2}+\mathcal{O}(a^{3})\,.$
(116)

For high spin values, a more accurate expression for these fluxes would be
found by expanding the fields to higher order in $a$. Doing so requires
boundary conditions to be imposed in order to determine the higher order
expansions of the fields and it not obvious how to do so (Tanabe & Nagataki,
2008). We can, however, find the expansion of $\dot{E}_{\rm BZ}$ to
$\mathcal{O}(a^{4})$ since the undetermined quantities in the higher order
expansion of the fields do not come into the energy flux at this order. We
find

$\displaystyle\dot{E}_{\rm BZ}$
$\displaystyle=8\pi\int_{0}^{\frac{\pi}{2}}\frac{1}{\sqrt{-g}}A_{\phi,\theta}^{2}\omega(\Omega_{H}-\omega)\sin^{2}\theta\,r_{\rm
H}\mathrm{d}\theta$
$\displaystyle=\frac{\pi}{24}C^{2}a^{2}+\bigg{(}\frac{56-3\pi^{2}}{1080}\bigg{)}\pi
a^{4}C^{2}+\mathcal{O}(a^{6})\,.$ (117)

These undetermined fields do, however, come into the second correction of the
angular momentum flux and so, without making further assumptions, we cannot
proceed further.

For the purposes of our sub-grid model, it is useful to express these fluxes
in terms of $\Omega_{\rm H}$ (Tchekhovskoy et al., 2010) and we do so using
the expansion

$a=4\Omega_{\rm H}-16\Omega_{\rm H}^{3}+64\Omega_{\rm H}^{5}+o(\Omega_{\rm
H}^{7})\,.$ (118)

Identifying $2\pi C=\Phi_{\rm BH}$ and reintroducing the factor of
$1/\sqrt{4\pi}$ that was absorbed into the Faraday tensor gives

$\displaystyle\dot{E}_{\rm BZ}=\frac{1}{24\pi^{2}}\Phi_{\rm BH}^{2}\Omega_{\rm
H}^{2}+\frac{(67-6\pi^{2})}{135\pi^{2}}\Phi_{\rm BH}^{2}\Omega_{\rm
H}^{4}+\mathcal{O}(\Omega_{\rm H}^{6})\,,$ (119) $\displaystyle\dot{J}_{\rm
BZ}=\frac{1}{12\pi^{2}}\Phi_{\rm BH}^{2}\Omega_{\rm H}+\mathcal{O}(\Omega_{\rm
H}^{3})\,.$ (120)

Thus far we have assumed a split-monopole configuration. A general magnetic
field configuration alters the fluxes only via the introduction of a factor
$\kappa$ which in the case of a split-monopole configuration is
$\kappa=1/(6\pi)$. In terms of $\kappa$ the fluxes are then

$\displaystyle\dot{E}_{\rm BZ}=\frac{\kappa}{4\pi}\bigg{(}\Phi_{\rm
BH}^{2}\Omega_{\rm H}^{2}+\frac{24(67-6\pi^{2})}{135}\Phi_{\rm
BH}^{2}\Omega_{\rm H}^{4}\bigg{)}+\mathcal{O}(\Omega_{\rm H}^{6})\,,$ (121)
$\displaystyle\dot{J}_{\rm BZ}=\frac{\kappa}{2\pi}\Phi_{\rm BH}^{2}\Omega_{\rm
H}+\mathcal{O}(\Omega_{\rm H}^{3})\,.$ (122)

Finally, we dimensionalise these equations and write them in terms of the
dimensionless magnetic flux

$\displaystyle\dot{E}_{\rm BZ}=\frac{\kappa}{4\pi}\phi^{2}_{\rm
BH}\dot{M}_{\rm BH,0}c^{2}\Bigg{[}$
$\displaystyle\bigg{(}\frac{a}{2(1+\sqrt{1-a^{2}})}\bigg{)}^{2}$
$\displaystyle+\frac{24(67-6\pi^{2})}{135}$
$\displaystyle\bigg{(}\frac{a}{2(1+\sqrt{1-a^{2}})}\bigg{)}^{4}\Bigg{]}+\mathcal{O}(\Omega_{\rm
H}^{6})\,,$ (123) $\dot{J}_{\rm BZ}=\frac{\kappa}{2\pi}\phi_{\rm
BH}^{2}\bigg{(}\frac{a}{2(1+\sqrt{1-a^{2}})}\bigg{)}\frac{GM_{\rm
BH}}{c}\dot{M}_{\rm BH,0}+\mathcal{O}(\Omega_{\rm H}^{3})\,,$ (124)

where we have used the fact that

$\Omega_{\rm H}=\frac{ac}{2r_{\rm g}(1+\sqrt{1-a^{2}})}\,.$ (125)