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

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.


INTRODUCTION
In the centres of many galaxy clusters the hot, X-ray emitting gas is observed to have rapid ( 1 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. 2019Russell et al. , 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. AGN jet feedback, however, is not confined just to these most ★ E-mail: rt421@cam.ac.uk (RYT) 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. 2019Mezcua et al. , 2018Chilingarian 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. 2012Mingo et al. , 2011Gallimore 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 1997Pringle , 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. (2020Liska et al. ( , 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.
There is a significant body of work that examines AGN jet launching on galaxy scales. These works investigate how the jet interacts with a multiphase interstellar medium (ISM) and how the gas structure in the vicinity of the black hole affects the resulting jet morphology on large scales. They also play a vital role in understanding the conditions under which jet launching may lead to enhanced star formation rates and those under which we would expect star formation to be quenched . To this end, some works consider the interaction of the jet with a single cold cloud (Antonuccio-Delogu & Silk 2008), others model a turbulent two-phase ISM (Asahina et al. 2017;Mukherjee et al. 2016;Wagner & Bicknell 2011) which may be confined to a disc (Cielo et al. 2018a;Mukherjee et al. 2018;Gaibler et al. 2012Gaibler et al. , 2011. Moving up to larger scales, much work has been done on jet propagation hydrostatic atmospheres and in idealised cluster setups, with the aim of better understanding the jet-ICM interaction (see e.g. Bambic 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;Sĳacki et al. 2015;Schaye et al. 2015;Crain et al. 2015;Dubois et al. 2012;Sĳacki 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. Due to resolution constraints, however, the jet energy must often be 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 & Sĳacki 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 & Sĳacki 2020;Bourne et al. 2019). To investigate these issues fully, however, self-regulated jet feedback models need to be considered (see e.g. Yang & Reynolds 2016b;Li & Bryan 2014b;Gaspari et al. 2011;Cattaneo & Teyssier 2007;Vernaleo & Reynolds 2006), but note that 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 -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 largescale 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 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.

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 steadystate, geometrically-thin -disc, as presented in Shakura & Sunyaev (1973). The black hole has a sub-resolution mass of BH and angular momentum J BH . Similarly, the disc is characterised solely by its mass, d , and total angular momentum, J d . The sub-grid disc-black hole system interacts gravitationally with the wider simulation via its dynamical mass, dyn = BH + d .
Our implementation of accretion via a thin, potentially warped -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.

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, , parameterised by , according to where s is the sound speed and is the disc scale-height. Whilst the value of this -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 material 1 from the inner edge of the -disc. If the rest-mass flux onto the black hole is BH,0 then the corresponding increase in black hole mass will be (acc) 1 Specifically the binding energy associated with the material.
where 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 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, d will change according to where in is the mass inflow rate from the surroundings onto the black hole-disc system and 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, BH,0 and the mass flux through the disc, . 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 where 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 ISCO in terms of the black hole spin) and j BH is the unit vector in the direction of the black hole angular momentum. Similarly, the disc angular momentum evolution is given by where L 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, , at the Eddington rate, Edd . Parameterising in terms a fraction Edd of Edd we can write = min( Edd , 1) Edd , where Salp is the Salpeter time and ES ≈ 0.4 cm 2 g −1 is the electron-scattering opacity. 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 where is the dimensionless spin parameter Theory constrains to lie in the range 0 ≤ ≤ 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 ≤ 0.998 (Thorne 1974). At low accretion rates we expect radiative cooling to become inefficient and the disc to thicken (Narayan & Yi 1994Abramowicz 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 Edd > (min) Edd = 10 −4 , as in Fiacconi et al. (2018).
Equation (3) admits the possibility that the accretion disc could grow without bound if in is consistently higher than . 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 For the derivation of the above expression for Fiacconi et al. (2018).

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 warp , is determined by a balance of the vertical warp propagation timescale, 2 , and the Lense-Thirring precession timescale, LT , where 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, J tot = J BH + J d is conserved, meaning that an equal and opposite torque must be felt by the black hole itself in order that J 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 where j BH and j 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 J tot and the second term leads to alignment with J tot (King et al. 2005). In general 1 and 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 1 and 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 1 and 2 as where GM is the gravitomagnetic timescale: (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 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 mass 2 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 j BH · j d > − BH /2 d (King et al. 2005), otherwise they will counter-align. Blandford & Znajek (1977) proposed a mechanism through which magnetic fields could extract rotational energy from the black hole.

The Blandford-Znajek mechanism
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 2 For the derivation of this warping mass see appendix A of Fiacconi et al. (2018).
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 where is the radial component of the magnetic field and d is the area element on the horizon then, in the limit of 1 and with fixed Φ BH , the rotational energy of the black hole can be extracted at a rate where g = BH / 2 is the black hole gravitational radius and is a parameter that depends only weakly on the magnetic field geometry. In this work we fix = 1/(6 ) 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 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 where the dimensionless magnetic flux is and ( ) is a dimensionless function of the black hole spin, given by This outward flux of energy leads to the evolution of black hole mass according to where, in analogy with the radiative efficiency, we have introduced a 'Blandford-Znajek efficiency' The corresponding outward angular momentum flux leads to evolution of the black hole angular momentum which we have defined in terms of an effective specific angular momentum The dimensionless magnetic flux that threads the black hole horizon, 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 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.

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 massloading factor, J (Bourne & Sĳacki 2017;Ostriker et al. 2010), a free parameter, which we define to be Throughout this work we choose J = 1 (Bourne & Sĳacki 2017).
In terms of the mass-loading factor, we can express the rest mass flux across the black hole horizon as The energy and momentum fluxes into the jet are where 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 It should be noted here that in our numerical scheme we inject the jet energy into a finite mass. Since this leads to numerical massloading, 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 subresolution velocity, 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.

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 and the disc mass The black hole angular momentum is given by and that of the disc is given by In terms of , the energy, mass and momentum fluxes into the jet are

NUMERICAL IMPLEMENTATION
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.

AREPO
We have incorporated this new black hole accretion and feedback prescription into the moving-mesh code (Springel 2010;Pakmor et al. 2016;Weinberger et al. 2020).
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.

Gas inflow properties
At the beginning of each black hole timestep we first estimate the gas inflow rate, in , and the specific angular momentum of this gas, L in , using a smooth particle hydrodynamic interpolation of the local mass flux onto black hole, calculated within the black hole smoothing length, ℎ BH 3 . 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 -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 where ∈ {1, 2, 3, 4} indexes the sector, and u are the density and velocity fields, respectively, and dS is the outward area element. F BH, is the numerical estimate of the mass flux evaluated at the position of the black hole for this sector and A is an effective area. We approximate the mass flux in each sector by summing over all gas cells that are radially inflowing and located within the relevant sector where is the mass density of the cell, , = (r − r BH ) · (u − u BH )/ r − r BH is the radial velocity of the cell and ( ) is a cubic spline kernel with compact support over ℎ BH and = r − r BH /ℎ BH .
We define the effective area to be A ≡ 2 in, where in, is a characteristic inflow radius calculated for each sector where is the volume of the gas cell. Similarly, the specific angular momentum of the inflowing gas in each sector is given by where L 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 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 in, = 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 L in using equation (45), but with the sum over all cells within ℎ 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 4 .
After this circularisation check has been carried out we are then able to combine the values of in, and L in, to get a global values for the mass inflow rate, in , and the specific angular momentum of inflowing material, L in , by summing over all sectors, with appropriate weightings. We store the values of in, 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 in and L 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 -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 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 in to ensure that the disc remains sub-Eddington by linearly interpolating equation (8) to ensure Edd < 1. If any of these checks lead to in being reduced then we also decrease the mass inflow rate in each sector (where in, ≠ 0) by the relevant amount.

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 4 We have verified that using 2-16 sectors in this calculation gives consistent estimates for the inflow properties. the predictor step and 'n+1', a value after the corrector step, this integration procedure corresponds to BH,temp

Black hole and disc angular momentum evolution
The angular momenta of the black hole and disc evolve according to equations (36) and (38). We implement this numerically by splitting the calculation into two steps. In the first step we evolve thedisc 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 where we are now using the subscript 'temp' to indicate the angular momenta after this step. The factor of J /(1 + 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, J /(1 + J ) Δ . Evolution due to the mass that flows onto the black hole from the inner edge of the disc, 1/(1 + J ) Δ , 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, J cons = J BH,temp + J d,temp is now conserved.
Having determined J 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 BH,n+1 = BH,temp + 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 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 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.

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 where indexes the properties associated with an individual gas cell and Δ in, = in, Δ .

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 -disc.

Jet injection
We implement jet injection using a procedure similar to that of the 'kinetic energy conserving jet' presented in Bourne & Sĳacki (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 Jet and height 2 ℎ 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). Jet varies such that the total mass within the jet cylinder, cyl , is constant 5 . Throughout this work we choose our target jet cylinder mass to be cyl = 0.5 M . ℎ Jet is determined by Jet such that the aspect ratio of the cylinder is fixed, giving a constant jet opening angle, Jet , that satisfies Following on from Bourne & Sĳacki (2017), we choose Jet /ℎ Jet = 3/2 which means that the cylinder has total volume (4 /3) 3 Jet . The cylinder is split into two halves, the north and the south (each with height ℎ 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 where is the initial mass of the cell, J is as determined by equation (41), the factor of a half accounts for the two halves of the jet cylinder and weight is the weighted sum of the masses of all gas cells in the relevant half-cylinder of the jet where is the (cylindrical) distance of the gas cell from the jet axis and is its height above the black hole equatorial plane. The kernel function has the form Momentum injection is carried out by giving gas cells a kick along the jet axis with magnitude where p and kin are the initial momentum and kinetic energy of the cell and Δ tot is the desired change in energy of the cell where J is as in equation (39). This single momentum kick, however, does not necessarily 5 We 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. ensure that the total energy of the cell changes by Δ 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, Δ therm , to ensure energy conservation, where 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).

Refinement
Our model makes considerable use of '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.

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 & Sĳacki (2015) which adaptively increases the mass and spatial resolution in the region surrounding the black hole.
Within a chosen refinement radius, ref , cells are forced to have radii 6 , , between two values, min ( ) and max ( ) which increase linearly with distance from the black hole, . Specifically, for a cell with ≤ ref , we ensure where Here, , max cell , and min cell are free parameters to be specified. It should be noted that our choices for these parameters differ from those suggested in Curtis & Sĳacki (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. 6 When we refer to the cell radius, we mean the radius of a sphere with volume equal to the cell volume:

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 & Sĳacki (2017). We briefly describe the relevant details of these schemes. Firstly, from the jet axis out to a (cylindrical) radius of Jet , the cell masses are forced below a spatially dependent maximum mass max cell ( ) which ensures that the jet cylinder is sufficiently well populated. max cell ( ) decreases as one gets closer to the jet axis according to where r is the (cylindrical) distance to the jet axis, , and are parameters to be chosen and is given by meaning that max cell Both this refinement scheme and the black hole refinement scheme detailed above only act within the 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 's mass refinement criterion (Springel 2010;Weinberger et al. 2020).
We identify jet material using a passive tracer which tracks the jet material, 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, , 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 ( ,J ≡ ,J / ) greater than 0.01 if their volume satisfies where max J 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 & Sĳacki 2017) and to limit the volume of cells as well as the gradients between neighbouring cells (Weinberger et al. 2017b).

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 where d / 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)).

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 ∼ 100 pc and masses 10 8−9 M (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.

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) where the bulge mass is * = 5 × 10 8 M and the scale radius is s = 100 pc. The gas disc is initialised with an exponential surface density profile where gas = 10 8 M is the disc mass and ℎ ≈ 70 pc is its scalelength. The initial vertical structure of the gas disc is determined such that it is in hydrostatic equilibrium and is initially isothermal with temperature 2 × 10 4 K.
We sample the stellar bulge with 2 × 10 6 collisionless particles of mass * = 250 M and with gravitational softening length * = 5 pc. The gas disc is initialised with 4 × 10 5 mesh generating points with target mass t = 250 M and gravitational softening length gas = 0.5 pc.
At the centre we place a black hole of mass BH = 10 6 M and gravitational softening length BH = 5 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 CGM = 10 −27 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. CGM = 10 −25 g cm −3 , which we refer to as the 'dense' CGM.
These initial conditions are then evolved for 200 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.

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, ,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 We use this tracer to asses the levels of mixing and entrainment of CGM material in the jet lobes.

Additional refinement
As discussed in Section 3.7, the mass refinement scheme in maintains gas cells at a chosen target mass, which we choose to be t = 250 M 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 & Sĳacki 2017).
Due to the density contrast between the circumnuclear disc and the CGM it is not necessarily the case that a global target mass t = 250 M 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 M . 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, piecewisedefined target mass that is specifically applied to the gas cells in the CGM 7 , targeting the increased resolution to a central cylindrical 7 Note 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. region (with the cylinder axis parallel to the -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, , to be CGM t

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, Edd 8 , 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 Edd -the 'low' power run, 10 −3 Edd -the 'medium' power run, and 10 −2 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 J = 10 5 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 & Sĳacki (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 scales 9 and, crucially, we are 8 Here Edd ≈ 1.26 × 10 44 erg s −1 refers to the Eddington luminosity of the initial black hole mass, BH = 10 6 M , 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. 9 In Bourne & Sĳacki (2017), the target jet cylinder mass is 10 4 M whereas in this work it is 0.5 M .
launching jets from a black hole embedded in a circumnuclear disc, whereas Bourne & Sĳacki (2017) launched jets from black holes embedded in a hydrostatic atmosphere.

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 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 where ram = 2 is the ram pressure and therm is the thermal pressure (Bourne & Sĳacki 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 / 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 -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 Myrs after the jet is launched. For the high power jet we find Mach numbers of ∼ 270 and ∼ 160 for the 'standard' and 'dense' jets, respectively, whereas those found in the 'low' power jets are significantly lower: ∼ 50 for the 'standard' run and ∼ 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 ∼ 150 pc in the 'low' power run to ∼ 400 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 & Sĳacki 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).

Comparison with analytic predictions
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)  uniform density CGM 10 : where h is the cross-sectional area of the bow shock at the end of the cocoon and h is the velocity of the head of the jet. Begelman & Cioffi (1989) assume that the thrust of the jet is given by Π J ∼ J / J where J is the jet power and J is the jet velocity. Bourne & Sĳacki (2017) used this model to validate the jets in their simulations, under the assumption that the jet thrust is given by Π J = J J /2 where J is the sub-grid jet mass flux, 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 is a constant which takes into account any momentum boost that occurs in the initial launching of the jet. Bourne & Sĳacki (2017) also expressed the cross-sectional area of the bow shock in terms of a working surface radius: h = 2 ws 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 be a constant is the assumption that the entire momentum flux at the base of the jet, J J /2, is what is felt at the head of the jet. In our analysis, we relax the assumption that is constant and instead of 10 Begelman & 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. assuming that ws scales linearly with time, we make a more general assumption that ws / 1/2 = + for some constants and . This is consistent with the fitting procedure in Bourne & Sĳacki (2017) but alters the interpretation of the result. In our case, it should be made clear that ws / 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) h ≈ Inserting ws / 1/2 = + then gives and integrating this then gives the jet length where 0 is the time at which the jet length is 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 J,0 = 500 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 pc or 2 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 and .
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 ws / 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 right 11 .
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 Myrs after the jet turns on (which corresponds to − 0 = 0.49 Myrs) after which its length stagnates. This behaviour is obviously not well fit by just a single pair of values, ( , ), 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 material 12 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 − 0 .
We have made use of '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 offaxis. Furthermore, the fact that the momentum flux at the head of the jet is time-variable justifies our interpretation that 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 ws / 1/2 , displayed in the bottom row of Fig. 4, should be linear.

Evolution of the jet lobes
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: where we choose the threshold jet fraction to be (thresh) J = 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 (thresh) J = 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 -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, respectively 13 .

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 13 Were 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. 10 2 LEdd Standard 10 3 LEdd 10 4 LEdd 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 -axis of each panel we have plotted blue ticks that are equally spaced in time.
kinetic and thermal components of this total lobe energy are shown by the orange and purple lines, respectively 14 .
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, the ratio of the lobe energy to the injected energy decreases with jet power, while in the 'dense' CGM runs there is not a clear trend.
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 pc this amounts to ∼ 0.32 for the 'high' power jet, ∼ 0.26 for the 'medium' power jet and ∼ 0.24 for the 'low' power jet. 14 We 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. 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 Myrs ( J = 874 pc) in the 'low' jet power case to 0.13 Myrs ( J = 462 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 CGM = 10 7 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 disc = 2 × 10 4 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 important 15 .
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, ∼ 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 ∼ 10 8 K. Whilst a non-negligible amount of disc material ends up in the cocoon in which we do find temperatures that exceed 10 8 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 15 It 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. the blue ticks on the -axis). Interestingly, we find that this ultimately leads to the energy in the lobe of this jet exceeding the total injected energy.

Lobe mass evolution
We now examine the evolution of the mass of the jet lobes in each of our simulations, shown 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.
Considering the runs in the 'standard' CGM case, it is apparent that the mass of these jets is entirely dominated by entrained disc material. 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, 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.
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.

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. We have 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.

Energy distribution
In Fig. 8, the total energy distribution is then further broken down into the kinetic energy, total thermal energy and thermal energy from entrained CGM material 16 . 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 and that the 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 shockheating 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 16 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.
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. In the 'low' power jet, however, the thermal energy from CGM material does becomes comparable to the total thermal energy of the jet lobes. 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. 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 z [pc] 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 (thresh) J = 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/ lobe ) d lobe,therm ( )/d at the recollimation shock then drops by a very large factor (of ∼ 18) 17 . 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) 18 . 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. In the case of the 'low' power jet, however, the thermal energy distribution 17 Note that in the 'standard' CGM case the results are much less sensitive to the adopted choice of (thresh) J as expected, given the much less developed turbulent cocoon. 18 One 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.
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 shockheating. 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. 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.

Mass distribution
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. 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. 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.

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 -disc is initialised with mass d = 10 3 M and angular momentum direction j d =ẑ, while the initial black hole mass and spin are BH = 10 6 M and = 0.2. This choice of -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 (36) and (39)), we evolve the simulations for 2 Myrs before turning on the jet to allow the sub-grid disc to come 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 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 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 -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 -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 -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. 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.
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 -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.

Evolution of the black hole properties
We begin our analysis by describing the evolution of the black holeaccretion 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 -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 -disc (right-hand panel).
In all runs the initial jet power is ∼ 10 40 erg s −1 , which is comparable to the 'low' power jets 19 which were presented in Section 5. Before the launching of the jet there is sustained mass inflow onto the -disc in all runs 20 . 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 ∼ 4 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 -disc. With no mass flowing onto the -disc to replenish the material accreted by the black hole or entrained by the jet, the mass of the -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 ∼ 1.7 Myrs and ∼ 4.0 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).
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--disc system. The amount of torquing is significant, with spin direction being changed by ∼ 30 • over the course of 10 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 -disc and the spindown 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 -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 -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 -discs are all initialised with a mass of 1000 M , but during the first 2 Myrs before the jet turns on, in all three runs the discs settle into quasi-equilibrium with mass of ∼ 400 M . Once the jets are launched, the -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 When considering the net angular momentum direction of the -disc we find that the discs in the 'vertical' jet cases remain directed along the -axis, as would be expected since the black hole spins remains approximately parallel to the -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 -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--disc system, albeit counteracted somewhat by sporadic inflow from the circumnuclear disc. It will be interesting to probe the effects of the initial -disc configuration in future work.

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).
In Fig. 11 we show slices in the -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 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 Myr the jet in the 'standard' CGM is already undergoing recollimation and at 10 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 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 ∼ 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 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 ∼ 5 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 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 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.

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 theplane (again, centred on the black hole), evolving in time with 3 Myr intervals from top to bottom. The main panel then shows a corresponding slice of the density field after a further 3 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 -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. J > 10 −3 ) and also has a -velocity greater than 50 kms −1 . After 10 Myrs we find that the outflowing material has mass ∼ 7 × 10 3 M 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 ∼ 3 and ∼ 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 -disc around ∼ 3 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 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 -disc the inclination angle decreases, with distinctive asymmetric features occurring in the outflow. After 10 Myrs, the jet has begun to emerge from the disc, as can clearly be seen in the inset temperature plot.

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 and cluster-scale 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) -disc.
Many jet studies use a simplified prescription for the jet power and direction and assume both to be fixed (e.g. Bourne & Sĳacki 2020;Bambic & Reynolds 2019;Ehlert et al. 2018;Bourne & Sĳacki 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. (2010Dubois et al. ( , 2012 adopt the Bondi-Hoyle-like accretion rate as a measure of the hot gas accretion rate, while (Yang & Reynolds 2016b;Li et al. 2017Li et al. , 2015Li & 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. 2018b). 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 & Sĳacki 2020Morsony 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 & Sĳacki 2015) and, further to this, Bourne & Sĳacki (2020 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 & Sĳacki 2020;Ehlert et al. 2018;Bourne & Sĳacki 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 & Sĳacki 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 & Sĳacki (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;Sĳacki 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 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. Previous studies considering the interactions between inclined AGN jets and the multiphase ISM of a galactic disc found evidence for morphologically similar outflows. Mukherjee et al. (2018) found that inclined jets launch sub-relativistic outflows along the galactic minor axis and similarly, Cielo et al. (2018a) found that the jet beams were deflected out of the disc plane. The fact that our simulations, and those mentioned above, find that AGN jets are able to drive outflows that closely resemble AGN winds is very intriguing and more work needs to be done to understand whether there are any observational signatures that could differentiate between these two driving mechanisms.

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 ab-initio 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 . Similarly, we do not consider special relativistic effects in our jet model as the velocities in all the jets considered in this work never exceed 0.2 . Any such effects will therefore have negligible impact on the propagation of our jets.
Our sub-grid disc is modelled as a thin -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 -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 disc 21 (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 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 21 Super-Eddington accretion is also likely to lead to the formation of a thick disc due to the larger optical depths (Abramowicz et al. 1988). 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;Sĳacki 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).

CONCLUSIONS
In this work we have developed a self-consistent sub-grid model for black hole accretion through a (warped) -disc and feedback in the form of a kinetic Blandford-Znajek jet within the moving mesh code . 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 selfregulated 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 -disc is drained. This leads to slower jet velocities, weaker shocks and colder material in the jet channel.
• The jet power is subject to self-regulation: as a consequence of the secular evolution of the circumnuclear disc, inflowing gas is eventually able to reach the -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 -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. , where Λ( ) = 3 + 2 ( ) ∓ √︁ (3 − 1 ( ))(3 + 1 ( ) + 2 2 ( )) .
The spin dependent radiative efficiency is then Finally, the specific angular momentum at the ISCO is 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 = BH = = 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 = /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 where is the stress-energy tensor, H = g (1 + √ 1 − 2 ) is the radial coordinate of the black hole horizon and is the determinant of , the metric tensor. Under the force-free condition reduces to the electromagnetic stress-energy tensor, EM , which is given by where is the Faraday tensor 22 which, in terms of the electromagnetic four-potential, , is Henceforth we will drop the subscript 'EM' in the electromagnetic stress-energy tensor.
The assumption of ideal MHD gives the constraint * = 0 , where * is the dual of the Faraday tensor. It then follows that , , = , , , and is a function of and therefore constant on magnetic surfaces (surfaces defined by = const). Thus, we define such that ( ) is a function that is also constant on magnetic surfaces. Imposing stationarity and axisymmetry then gives = = 0. The remaining non-zero components of the Faraday tensor can then be written in terms of the three variables , and as and we also identify 22 We have absorbed a factor of 1 √ 4 into the definition of the Faraday tensor.
The relevant components of the stress-energy tensor, evaluated on the black hole horizon, are then The evaluation of equations (B1) and (B2) therefore require the determination of the fields , and . This is done by solving the equation of motion EM ; = 0 .
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 . For this = 0 case, the equation of motion reduces to where we use the superscript (0) to indicate solutions to this nonspinning case and L is the linear operator L ≡ 1 sin 1 − 2 + 1 2 1 sin . (B14) The solutions for a split-monopole field geometry are where is a constant. Perturbing this to the next highest order in , the fields take the form ( , ) = − cos + (2) ( , ) 2 + O ( 4 ) , is even in since the radial magnetic field geometry should not depend on the direction of the rotation. Symmetry arguments further constrain and to be odd functions of .
To O ( 2 ) the equation of motion is where ( , ) is a source function that depends on (1) and (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 (1) and (1) , and can thus be found by inverting equation (B21) using the known Green's function for the linear operator. Altogether this gives ( , ) = − cos + ( ) cos sin 2 2 + O ( 4 ) , ( , ) = 1 8