Non-isotropic feedback from accreting spinning black holes

Active galactic nuclei (AGN) are massive black holes (BHs) caught in the act of accreting gas at the centre of their host galaxies. Part of the accreting mass is converted to energy and released into the surrounding medium, in a process loosely referred to as AGN feedback. Most numerical simulations include AGN feedback as a sub-grid model, wherein energy or momentum (or both) is coupled to the nearby gas. In this work, we implement a new momentum-driven model in the hydrodynamics code GIZMO, in which accretion from large scales is mediated by a sub-grid accretion disc model, and gas particles are stochastically kicked over a bi-conical region, to mimic observed kinetic winds. The feedback cone's axis can be set either parallel to the angular momentum of the gas surrounding the BH or to the BH spin direction, which is self-consistently evolved within the accretion-disc model. Using a circumnuclear disc (CND) as a test bed, we find that (i) the conical shape of the outflow is always visible and is weakly dependent on the launching orientation and aperture, resulting in comparable mass inflows and outflows; (ii) the cone's orientation is also similar amongst our tests, and it is not always the same as the initial value, due to the interaction with the CND playing a crucial role in shaping the outflow; and (iii) the velocity of the outflow, instead, differs and strongly depends on the interplay with the CND.


INTRODUCTION
Supermassive black holes (BHs) with masses in the range 10 6 -10 10 M are typically found at the centre of massive galaxies (Kormendy & Richstone 1995). They represent the most energetic and compact sources in the Universe, and are typically observed as active galactic nuclei (AGNs), through the radiation emitted by the material accreting on to them. According to the commonly accepted picture, BHs form in the early Universe and then grow with cosmic time via accretion and mergers, releasing an enormous amount of energy in the process (dubbed 'AGN feedback'), that can significantly affect their surroundings. Because of their powerful feedback, AGNs are usually advocated to explain the quenching of star formation in the most massive galaxies (Silk & Rees 1998;Di Matteo et al. 2005), and have been also proposed as potential candidates to solve other open problems, such as the source of hydrogen reionization  2009), the too-big-to-fail problem (Garrison-Kimmel et al. 2013), and the removal of gas from massive disc progenitors (e.g. Peirani et al. 2012), although a clear consensus has not been reached yet.
AGN feedback, because of its role in the BH mass growth self-regulation, is also invoked as a mechanism to explain the emergence of the observed correlations between the central BH and the velocity dispersion of the central stellar spheroid (e.g. Ferrarese & Merritt 2000;Gebhardt et al. 2000;Häring & Rix 2004;Gültekin et al. 2009;Kormendy & Ho 2013) or its stellar mass (e.g. Magorrian et al. 1998;McConnell & Ma 2013;Reines & Volonteri 2015). In particular, different scenarios have been proposed to date to explain these correlations, wherein (i) the galaxy grows first and the BH catches up later in a supply-limited fashion, (ii) the BH is initially overmassive and the stellar component is built-up only after the AGN feedback-regulated regime has been reached, or (iii) stars and BH co-evolve along the correlation in a self-regulated fashion (see reviews by e.g. Volonteri 2012; Kormendy & Ho 2013). Despite the increas-ing amount of data available, and the information obtained at high redshift, a clear consensus about how and when these correlations arose is still missing.
In this respect, AGN feedback can play a crucial role, and the observation of jets and large-scale winds represents an extremely important evidence of this mechanism at play. In particular, fast nuclear (sub-pc) outflows (Tombesi et al. 2013;Nardini et al. 2015) as well as galaxy-scale (kpc) winds (Feruglio et al. 2010;Rupke & Veilleux 2011;Sturm et al. 2011;Greene et al. 2012;Maiolino et al. 2012;Liu et al. 2013;Cicone et al. 2014;Harrison et al. 2014) allow one to infer how much energy and momentum are injected in the BH surroundings and transported up to kpc scales, affecting the evolution of the system, and can be considered as the indirect evidence of the existence of such feedback processes.
In the X-ray band (Gofford et al. 2015), ultrafast outflows (with velocities larger than 10 per cent of the speed of light) are observed in 40 per cent of the bright local population of AGNs. Deep absorption lines associated with highly ionized iron (Reeves et al. 2003;Tombesi et al. 2013;Gofford et al. 2015) are also typically found to be blueshifted with velocities in the range ∼10 3 -10 4 km s −1 , suggesting the presence of outflowing material directly connected to the central BH energy injection on sub-pc scales. From the detection fraction of 0.4, an opening angle can also be derived for these fast pc-scales winds (as done in e.g. Gofford et al. 2015), corresponding to a region encompassed by the outflow which subtends almost half of a sphere.
Recently, massive galaxy-size cold molecular outflows (up to 10 10 M in molecular hydrogen) have been detected with velocities in the range ∼10 2 -10 3 km s −1 at far-infrared and millimetre wavelengths (e.g. Feruglio et al. 2010;Fischer et al. 2010;Cicone et al. 2014Cicone et al. , 2018García-Burillo et al. 2014;Sirressi et al. 2019). In addition, there are examples of winds associated with molecular gas, with speeds of the order of 10 3 km s −1 , showing both collimated and un-collimated geometries at hundreds of pc scales (see e.g. Combes 2017, and references therein).
Theoretically, several works have addressed the role of BH feedback in the emergence of galaxy-BH correlations and in galaxy quenching (e.g. Silk & Rees 1998;Fabian 1999;King 2003), showing that, when the luminosity produced by the central object couples with the nearby matter, a fast nuclear wind is generated, which can interact with the interstellar medium (ISM) of the host, sweeping up material and driving large-scale outflows able to suppress star formation. In addition, by accounting for the energetics of the shock-generated shell after an AGN feedback event, the scaling relation between the stellar spheroid's velocity dispersion σ and BH mass MBH can be derived (King 2003) either in the fully energy-conserving regime (MBH ∝ σ 5 ) or in the momentum-conserving case (MBH ∝ σ 4 ), bracketing extremely well the observed correlation, i.e. MBH ∝ σ 4.38 (Kormendy & Ho 2013).
However, given the complexity of the AGN feedback problem and its interplay with the galaxy host, involving a non-linear coupling of multiple physical processes on extremely different scales (from the accretion disc up to the cosmological environment), a careful combination of observations, analytical modelling, and numerical simulations is needed. Over the last decade, many numerical hydrodynamic simulations have tackled this problem, focussing on different scales. For instance, at small scales, several studies have investigated the effects of radiation on the accretion itself and the wind/jet driving with radiation-hydrodynamics -in some cases, also including magnetic fields and general relativistic effects (e.g. Park & Ricotti 2012;Sądowski 2016;Jiang et al. 2019). At larger scales (and lower resolution), where the accretion disc scales cannot be properly resolved, effective sub-grid models were developed, mostly to address the effects of AGN feedback on the galaxy ISM and the BH dynamics (see e.g. Hopkins et al. 2007;Choi et al. 2012;Capelo et al. 2015;Zubovas et al. 2016;Anglés-Alcázar et al. 2017;Souza Lima et al. 2017;Lupi et al. 2019). Finally, at cosmological scales, where the spatial resolution is too low to properly resolve the galaxy ISM, AGN feedback is mainly modelled via an effective coupling efficiency of the accretionpowered luminosity with the gas, calibrated to reproduce the observed local correlations, as in e.g. Schaye et al. (2015); Dubois et al. (2014) and Vogelsberger et al. (2014).
Although most of the current state-of-the-art cosmological simulations have been shown to correctly reproduce many galaxy and quasar properties, we are still far from a proper understanding of how the AGN feedback actually couples on different scales and at different resolutions. For example, alternative explanatory scenarios for the correlations we observe are possible: Anglés-Alcázar, Özel & Davé (2013) and Anglés-Alcázar et al. (2015), using the sub-resolution accretion rate estimator based on gravitational torques by Hopkins & Quataert (2011), showed that feedback-regulated growth is not required for the emergence of scaling relations between the BH and its host galaxy.
In this work, we build upon the BH accretion and feedback models by Anglés-Alcázar et al. (2017) and introduce a self-consistent model of a sub-resolution Shakura & Sunyaev (1973) accretion disc that evolves via accretion from the surroundings (on resolved scales), draining due to the gas accretion on to the central BH, and wind production (i.e. outflows). We implement our model in the publicly available code gizmo 1 (Hopkins 2015), descendant of gadget2 and gadget3 (Springel 2005;Springel et al. 2008), and validate it via a suite of circumnuclear disc (CND) simulations. We also explore the role of the central BH spin, by coupling our sub-grid model with the self-consistent spin evolution model presented in Cenci et al. (2021).
The paper is organized as follows: in Section 2, we describe our accretion disc model and the coupling with the hydrodynamic code; in Section 3, we validate the model in a controlled environment; and in Section 4, we draw our conclusions.

IMPLEMENTATION
In this section, we describe the sub-grid model for BH accretion and feedback employed in this work, which builds upon the one presented in Anglés-Alcázar et al. (2017), where we add the sub-resolution modelling of an accretion disc and allow the outflow to have different shapes and directions. In the following, we provide the general description of the model, whereas all specific values of the parameters are given in Section 3.

Accretion from large scales
Since galaxy-and cosmological-scale simulations are unable to properly resolve the accretion disc around BHs (of the order of the Schwarzschild radius), accretion is usually modelled via semi-analytical prescriptions that link the resolved inflow rate to the actual accretion on to the BH.
In gizmo, the mass transfer rate from resolved scales onto a BH particle (corresponding to a sub-resolution 'BH+accretion disc' system) can be estimated using a variety of different prescriptions. In this work, we assume that the mass transfer rate on to the unresolved system is determined by the Bondi-Hoyle-Lyttleton formula (hereafter Bondi;Hoyle & Lyttleton 1939;Bondi & Hoyle 1944;Bondi 1952) as implemented by Springel et al. (2005) where ρ and cs are the density and the sound speed of the gas, respectively, evaluated as mass-weighted average quantities using the BH neighbouring gas particles, v is the velocity of the BH relative to the centre-of-mass velocity of the gas within the BH kernel length, MBH = 10 6 MBH,6 M is the BH mass, G is the gravitational constant, and αacc is a dimensionless parameter usually employed in ∼kpcresolution simulations to boost the accretion rate in the aim at correcting for the unresolved dense gas in the ISM (Di Booth & Schaye 2009). Due to the resolution of our tests, αacc is set to 1. The local properties of the gas surrounding the BH are averaged via a smoothing kernel defined as the region encompassing an effective number of neighbours N ngb,BH = fBHN ngb , where N ngb is the user-defined number of neighbours used for hydrodynamics and fBH ≥ 1 is a free parameter (usually chosen to ensure a good covering of the region surrounding the BH). In order to prevent the BH from unphysically accreting gas at large distances, the code allows to define a maximum kernel size used for accretion, i.e. the maximum accretion radius raccr,max, by means of an additional free parameter.

The original model
In the model by Anglés-Alcázar et al. (2017), two distinct masses are associated with the BH particles, dubbed the 'physical' and the 'dynamical' mass, respectively. The former is evolved by any chosen sub-grid accretion model in order to follow the BH mass growth, whereas the latter is the mass that is used to calculate the BH's gravitational force and dynamics.
The BH physical mass, MBH, evolves at each time-step as wherė Macc−BH is the accretion rate on to the BH (see below), ∆t is the BH time-step in the simulation, and η = 0.1η0.1 is the radiative efficiency, which can be kept constant during the simulation by setting its value in the parameters file of the code, or let free to vary self-consistently with the properties of the sub-grid accretion disc, as done in, e.g. Cenci et al. (2021). In the original model, the accretion on to the BH is not mediated by an accretion disc. Therefore, the accretion rate,Ṁacc−BH, is assumed to be a fixed fraction f of the large-scale inflow rate:

Mass conservation imposeṡ
Min =Ṁacc−BH +Ṁout, whereṀout is the outflow rate. Hence In a radiatively driven-wind scenario, gas is accelerated by radiation pressure, with the photons being produced during the accretion process yielding a total (bolometric) luminosity L bol = ηṀacc−BHc 2 , where c is the speed of light in vacuum. The outflow rate is then given by the definition of outflow momentum flux in units of L bol /c (hereafter "momentum-loading"): where vout is the outflow velocity's magnitude. The quantities vout and p b are free parameters of the model and can be set in the parameters file of the code. Equations (6) and (7) provide a definition for f : Once the physical BH mass growth has been computed, a set of gas particles are stochastically selected within the BH kernel, assigning to each neighbour a probability where MBH dyn is the BH dynamical mass at the beginning of the time-step, wj is the kernel weight, m k is the gas particle mass, and the sum is over the particles that have already been selected in the same time-step. In the case M BH > MBH dyn , the factor 1/f accounts for the mass that has to be ejected. After being selected, a gas particle has a fraction f of the mass removed and transferred to the BH dynamical mass, such that, stochastically, Hence, the BH dynamical mass keeps track of the discrete accretion, i.e. it increases by an amount equal to each accreted particle in the simulation (or a fraction of it). The physical mass, on the other hand, keeps track of the quasicontinuous accretion (that evolves with the BH time-step of the simulation, ∆t).
In the case M BH < MBH dyn , we only change the momentum of the selected particles, and the BH's dynamical mass does not change. This ensures mass conservation and, at the same time, that the dynamical mass follows the physical mass (on average).
In both cases, a velocity with magnitude vout is added to the velocity of each particle that has been selected in a specific direction, hence giving it a "kick" in that direction. Concerning the specific prescription for the direction, the implementation by Anglés-Alcázar et al. (2017) assumes that the kicks can be given either along the radial direction connecting the gas particle to the BH ('radial outflow') or along the angular momentum vector of the gas particles within the BH kernel length ('collimated outflow').
We note that, in this model, the overall amount of ejected mass per time-step is stochastically equal toṀout∆t, andṀout is defined once f -which depends only on free parameters and remains constant -has been provided.

The sub-grid model
In our implementation, we add a self-consistent subresolution model of an accretion disc. Since mass must be conserved for both resolved and unresolved structures, we write the rate of change of the disc mass, M disc = 10 4 M disc,4 M , aṡ where the inflow rate is computed from resolved large-scale quantities according to Equation (1), and the outflow rate is computed as in Equation (7). The accretion rate on to a BH is computed aṡ Eddington (1916) accretion rate, mp the proton mass, and σT the Thomson cross section. The Eddington ratio f Edd is given by (Fiacconi et al. 2018) where aBH = cJBH/(GM 2 BH ) is the BH spin parameter, J disc is the angular momentum magnitude of the accretion disc, α = 0.1α0.1 is the disc viscosity parameter, and JBH is the magnitude of the BH angular momentum, defined as JBH = JBH ·jBH, where jBH is its direction. Furthermore, note that, in our model η is computed as a function of aBH. We further impose that f Edd ≤ 1 and, at every time-step, we limitṀin to prevent the sub-grid disc from becoming self-gravitating (for more details on how α, η, and the limitation onṀin are computed, see Cenci et al. 2021). Finally, at any given time-step, if the mass subtracted from the disc is greater than the disc mass itself, the disc mass is set to zero. 3

Stochastic accretion and feedback
We evolve the BH physical mass as in Equations (2) and (3). We additionally follow the accretion disc mass that is updated at each time-step as whereṀ disc is given by Equation (12).
The stochastic selection probability in Equation (9) must therefore be modified as In the case M BH + M disc > MBH dyn ,Ṁout∆t accounts for the mass left to the particles after accretion, that will constitute the outflow mass. The mass fraction that has to be removed from each particle is consistently re-defined as where N is the total number of selected particles. The code removes this mass fraction from each particle and adds it to MBH dyn , as in Equation (10), such that, stochastically, Next, each particle is kicked with its remaining mass, which is then just ∆Mout/N . We stress that our redefinition of the accreted fraction from each particle is variable on a time-step basis, because it depends on the number of selected particles, which in turn depends on both the amount of mass that has to be accreted and the amount launched in the outflow. This redefinition is necessary because of the presence of the sub-grid accretion disc. When an accretion disc is present, the constant f in the previous model is not valid any longer, sinceṀout is defined by the accretion rate on to the BH, and not by that on to the disc. Hence, f is defined only onceṀout has been computed starting froṁ Macc−BH. However, the particles are selected such that the whole amount of ejected mass per time-step is stochastically equal toṀout∆t.
In the case M BH + M disc < MBH dyn , ∆MBH dyn = 0 and we only change the momentum of the selected particles.

Conical feedback
In order to better reflect observations of radiatively driven winds (see e.g. Gofford et al. 2015;Combes 2017), we implement an outflow with a biconical shape, whose axis can be either fixed in time (to a given arbitrary direction) or evolved during the simulation (e.g. following the angular momentum of the gas surrounding the BH). We further model the case in which the cone's axis is set parallel to the evolving BH spin, given by the prescription described in Cenci et al. (2021).
For each particle within the BH kernel that receives a kick, we randomly sample the kick direction within a cone of semi-apertureθ aligned with a chosen direction (either fixed or time-dependent). Hereafter, we will denote this direction, along which feedback is exerted, using a spherical coordinates system (θ, φ), respectively the polar and azimuthal angle. Note that, in our model, all the particles within the BH kernel are eligible to be kicked: only the direction of the kick is extracted such that the angles distribution describe a cone in the velocities space. This translates into a conical outflow assuming the shape of a cone (in the positions space) with the base corresponding to the BH kernel. This approach has the advantage of always finding particles to kick, unlike other schemes where, e.g. only particles within the cone are selected, possibly resulting in missing feedback (Barai et al. 2016).
In order to produce a symmetric outflow with respect to the plane perpendicular to the chosen axis, the algorithm then proceeds in two steps. First, it computes the radial position of each particle with respect to the BH, r, and the scalar product of this vector and the randomly extracted velocity vector, r · v. If the latter quantity is positive, then Table 1. Resolution of our simulation suite. For each particle type, N is the number of particles, m (or M ) the particle mass (in M ), and ε the gravitational softening (in pc). While the stellar quantities are fixed during the simulations, the gas and BH particle masses (as well as the number of gas particles) evolve, as a result of accretion events, hence we report here the initial values (with the subscript 0; also, M BH,0 + M disc,0 = M BH dyn ,0 , or M BH,0 = M BH dyn ,0 when the accretion disc model is not used). The softening is fixed for stars and the BH, whereas fully adaptive softening is employed for gas, whose minimum is here reported.

TESTS
In order to study the accretion and feedback processes in detail, we chose a test setup that allowed us to have a wellcontrolled environment, with the minimum amount of active sub-grid physics. The tests were performed employing the initial condition prescriptions from Lupi et al. (2015). The setup has three components: • a stellar spherical structure described by a Hernquist (1990) profile with total mass M b = 5 × 10 8 M and scale radius r b = 100 pc, defined by the radial density profile where r is the spherical radial coordinate.
• a gaseous CND with MCND = 10 8 M and scale radius RCND = 50 pc, defined by the surface density profile where R is the cylindrical radius. The vertical component and velocity field were computed as described in Lupi et al. (2015), using an iterative procedure that takes into account the global potential of the three components (stars, gas, and BH) and ensures vertical hydrostatic equilibrium. 4 In our tests, the CND initially lays in the xy-plane and, during the simulations, never significantly changes its orientation. The direction (θ, φ)=(0, 0) identifies the positive z -axis, whereas (θ, φ)=(π/2, 0) identifies the positive x -axis. Table 2. Summary of the parameters adopted in our simulations. The subscript 0 refers to quantities evaluated at initialization, i.e. at time t = 0. Moreover, the following parameters are the same for all tests: p b = 1, η 0 = 0.1 and, for the accretion disc model in the last four runs, a BH,0 = 0.5, M disc,0 = 5 × 10 4 M , f Edd,0 = 5 × 10 −3 , and R circ /Rsg = 0.5. • a central BH with initial mass MBH,0 = 10 7 M .

Simulation nameθ
Our tests were conducted at two different resolutions, as detailed in Table 1. We compared the high-and lowresolution results and verified that there is no significant difference in the system's evolution. Hence, in the following, we present only the high-resolution results.
The gas was initialized following an ideal-gas equation of state with γ = 5/3 and a uniform temperature of 2 × 10 4 K. Since, as described in Lupi et al. (2015), some spiral arms have to be dissipated and the disc undergoes an initial re-adjustment, we evolved our model in isolation with only gravity and hydrodynamics active for 20 Myr and reset the time to 0 at the end of this relaxation period. Afterwards, the tests with our model were performed with γ = 7/5, to effectively mimic a mild cooling and favour gas accretion towards the centre, without actually using a cooling model (see e.g. Dotti et al. 2009).
In all our simulations, we set N ngb = 32 and fBH = 3. We further assumed no boost factor in the Bondi accretion rate (i.e. αacc = 1 in Equation 1) and fixed raccr,max = 10 pc, which, in our setup, is ∼10∆x, where ∆x is the mean interparticle distance amongst gas particles.
The remaining parameters of our model vary depending on the test and are therefore quoted in the following sections.

Comparing different feedback geometries
Before testing our new accretion model (described in Section 2.3 and in Cenci et al. 2021), we validated our new feedback model (presented in Section 2.4), using the original accretion scheme by Anglés-Alcázar et al. (2017) and comparing their radial and collimated feedback geometries to our conical feedback prescription, for two different values of the kick velocity's magnitude: vout = 10 3 and 10 4 km s −1 (see Table 2 for a summary of our tests). In this case,Ṁacc−BH is computed using Equation (4), η is fixed to 0.1 in the parameters file, and the feedback reference direction (for the collimated and conical cases) is computed on a time-step basis as parallel to the angular momentum of the gas within the BH kernel (hereafter Jgas = Jgas · jgas, where Jgas and jgas are the magnitude and direction of the vector, respectively). Moreover, we assume a momentum-conserving scenario, by setting p b = 1. Hence, the values of f as defined is Equation (8) are f = 0.032 and 0.25 for vout = 10 3 and 10 4 km s −1 , respectively. The cone semi-apertureθ is set to π/4. Figure 1. Since the BH mass changes by less than one per cent over the whole 10 Myr of the simulations, the change in the Eddington ratio is approximately proportional to the BH accretion rate. At the very beginning of the simulations, accretion is quite high, because the BH is surrounded by a dense gas region. As soon as feedback becomes efficient, the accretion rate quickly drops. After this initial transient, which lasts ∼0.05 Myr, we have two different scenarios, depending on the kick velocity's magnitude.

The Eddington ratio for these tests is shown in
When vout = 10 3 km s −1 , for the first ∼4 Myr the collimated feedback has the higher accretion rate, that in the conical case is slightly lower, and that in the radial case is much lower. This is due to the fact that the accretion increases as the interaction with the CND decreases: when particles are mostly pushed towards the CND (as in the radial model), accretion is impeded and therefore is lower than when particles are mostly pushed in the vertical direction (as in the collimated and conical cases), leaving the accretion flow from the CND almost unaltered. Furthermore, the accretion rate in the conical case is slightly lower than in the collimated case because as the semi-aperture increases there are more particles interacting with the CND. After around 4 Myr, the pattern becomes more complex. This is due to the formation of asymmetric gas overdensities (visible in the face-on density panels of Figure 2, shown at 5 Myr) and the BH starts wandering as a result of the gravitational interactions with such overdensities. At later times, the BH can get close enough to these denser regions, such that the accretion rate and, consequently, feedback increase. This in turn pushes away the gas and decreases the accretion rate. This behaviour repeats, but it is heavily affected by where and when the overdensities form, which in turn is also influenced by the kicks that are intrinsically stochastic.
When vout = 10 4 km s −1 , the fraction f defined by Equation (8) is almost one order of magnitude higher than in the previous case (for a fixed p b and η). Hence, to have the same accreted mass, fewer particles are needed. As a result, the vout = 10 4 km s −1 case is more stochastically affected by fewer, higher-velocity particles that periodically . Note the different y-axis range in the two panels. In the case with vout = 10 3 km s −1 , in the first ∼4 Myr the Eddington ratio for the conical and collimated case is higher than in the radial case, where the larger interaction with the CND more effectively suppresses the accretion. In the following ∼6 Myr (and during the entire simulation with vout = 10 4 km s −1 ), the main effect is that the accretion increases and decreases repeatedly due to the interaction between BH and overdensities that form during the simulations, also as a result of the stochastic kicks.
lower the accretion by a significant amount, whereas in the other case (vout = 10 3 km s −1 ) the pattern is smoother. Each single particle is able to heavily modify the gas distribution around the BH in a specific direction in the CND plane. All the three tests with different kick's velocity magnitude exhibit the behaviour in which the accretion rate increases and decreases repeatedly due to the interaction between BH and overdensities. We note that the accretion rate in the vout = 10 4 km s −1 case is typically lower than that in the vout = 10 3 km s −1 case. It is however hard to generalize this result, given the short time-scale and the extremely idealized setup of the simulation, which result in a negligible BH mass growth. We plan to investigate such trend in future works. Figure 1 shows thatṀacc−BH/Ṁ Edd 0.01 even if we are using a feedback model that aims to represent a fast wind in a radiatively efficient accretion phase. However, the accretion rate depends mostly on the properties of the CND, which have been chosen to be in quasi-equilibrium, in order to prevent disc fragmentation, the need for additional physics such as star formation and its associated feedback. This choice allowed us to isolate the effect of the AGN feedback prescription, at the cost of having low accretion rates.
We verified that the kicked particles, for the collimated and conical cases, define a conical shape far from the launch region for the entire simulation, as illustrated in the bottom panels of Figure 2 for the vout = 10 3 km s −1 case, where the particles that have been given a kick are highlighted with the large dots. The semi-aperture parameter determines the large-scale particle direction dispersion, because it reflects the small-scale launch directions distribution. However, as soon as the wind leaves the high-density region of the CND, where gas pressure helps keeping it collimated, the gas expands in the radial direction, forming a cone (also in the collimated case).
Two more quantities must be provided to the accretion disc model to be initialized properly: we set the initial Eddington ratio 5 f Edd,0 = 5 × 10 −3 and the ratio of the circularization to self-gravitating radius Rcirc/Rsg = 0.5 (for more details, see Cenci et al. 2021).
Hence, the coordinates of the initial BH spin direction  2017), and BH feedback with vout = 10 3 km s −1 and radial (left-hand column), collimated (central), and conical (right-hand) geometry. First row: gas surface density (viewed face-on); the red, blue, and yellow dots indicate the position of the BH, the centre of mass of the gas, and the centre of mass of gas+stars, respectively. Second row: gas particles' positions map (viewed face-on); the particles that received a kick are highlighted with the large dots, whose colour code reflects the particle's z coordinate (greater values are lighter); the red circle shows the extent of the BH accretion length, whereas the red arrow shows the projection of the angular momentum direction of the gas in the BH kernel, arbitrarily normalized to be clearly visible. Third and fourth rows: same as the first and second rows, respectively, but with the CND seen edge-on (omitting the centre of mass positions) and with the gas velocity field overlaid to density map. In the radial case, the particles are heavily scattered when kicked along the disc plane, and a clear preferential direction of the particles after their launch is not recognizable. On the other hand, in the collimated and conical cases, after the particles have been launched, they define and maintain a conical shape. When the conical prescription is adopted, we can control the scatter in the resulting directions. The conical shape is still visible in the collimated case because, even if the particles are always launched in the perpendicular direction to the disc plane, without any direction distribution, the gas can freely expand in the low-density region above/below the disc, where there is not enough pressure to keep the wind confined.
For simplicity, we start by exploring how the choice of the cone's axis direction and semi-aperture affects the resulting outflow in a series of tests in which η andṀacc−BH are taken as an input from the accretion disc module, but the cone axis direction is fixed in time (although the BH spin evolves). We are particularly interested in assessing how different relative orientations between the feedback reference Figure 3. Same as Figure 2, but for the tests described in Sections 3.2 and 3.3. From left to right: along_z_45, along_jBH0_45, along_jBH0_10, full_model. Considering along_z_45, the kicked particles at large scales (large dots in the bottom panels) define a conical shape, which is maintained for the whole duration of the tests. Tests along_jBH0_45 and along_jBH0_10 show no significant dependence of the result on the semi-aperture. Regarding full_model, the red arrow highlights that the BH spin, at ∼5 Myr, is directed along the disc plane, hence the large-scale outflow shape is more affected by the interaction of the wind with the disc (the particle distribution looks more cluttered, especially near the centre).
direction and the CND symmetry axis affect the outflow properties.
(ii)θ = π/4 and cone's axis direction constant throughout the test and equal to the value of jBH0 given by Equation (20).
(iii)θ = π/18 and cone's axis direction as in the previous case.
We forced the direction to be a definite value dur-ing these simulations, which is a reference-frame-dependent choice, because the CND starts in the xy-plane and we verified that its orientation does not change significantly during the simulation under the effect of feedback. In all three cases, as well as in the "full model" test described in the following section, we assume a momentum-conserving scenario (i.e. p b = 1) and a kick velocity magnitude vout = 10 3 km s −1 .
The test along_z_45, in which feedback is exerted along the z-axis, is similar in nature to the conical-feedback case of the previous section, except for the mass rates, which are now computed using our sub-grid accretion disc model. The first column of Figure 3 looks very similar in every aspect to . From top to bottom, respectively: along_z_45, along_jBH0_45, along_jBH0_10, and full_model at 5 Myr. Left-hand panels (polar plots):θ (blue line) and instantaneous feedback direction (purple line) as a function of radius, plotted for the snapshot. Central panels: xy-projections of the upwards (left) and downwards (right) outflow particles positions. Right-hand panels: xz-projections of the upwards (left) and downwards (right) outflow particles positions. The colour code reflects the particle z coordinate (lighter are further from the xy-plane). The polar plots show that the direction of the resulting outflowθ varies with the distance from the BH and is not exactly equal to its launch value (shown as the purple line). The visual appearance of the outflow (cartesian panels) is also quite similar.
the last column of Figure 2, as expected, since the accretion model (the only difference between the two cases) affects the mass rates (hence the number of selected particles) but it has no relation with the outflow geometrical features. Moreover, the kicked particles at large scales (large dots in the bottom panels) define a conical shape, which is maintained for the whole duration of the tests.
We further defined the outflow as the subset of gas particles which have v radial > 2× v tangential (where v tangential is the mean tangential velocity over the CND and maintains an approximately constant value of ∼85 km s −1 ), and that are at a radius greater than the BH accretion length. We divided the outflow in spherical shells centred on the BH position and, for the particles belonging to each shell, we computed the median position, which corresponds to a direction (θ) in spherical polar coordinates. This direction, computed shell by shell, is therefore dependent on the distance from the BH.
Considering the blue line in the polar plot of the top row of Figure 4, the outflow directionθ deviates at every radius from being exactly on the z -axis, because of two concurrent effects: the particle scattering effects when interacting with the surrounding disc and the pressure gradients that develops due to hydrodynamics, since the wind is expanding in the low-density region above/below the disc.
Comparing tests along_jBH0_45 and along_jBH0_10 (see the two central columns of Figure 3 and the two central rows of Figure 4), which only differ with each other in the value of the semi-apertureθ, we note that this difference has no significant effect on the resulting outflow. As in the case of test along_z_45, the direction of the resulting outflowθ varies with the distance from the BH and is not exactly equal to its launch value (shown as the purple line in Figure 4).
In Figure 5, we show the velocities' magnitude of the outflowing gas. In test along_z_45 (top-left panel), the speed reaches more than ∼550 km s −1 , whereas in tests along_jBH0_45 and along_jBH0_10 (respectively, top-right and bottom-left), it reaches ∼400 km s −1 . This effect can be explained considering that the kicked gas particles shock 3. An inflow from the disc plane is quite evident, but note that its velocity is two orders of magnitudes lower than that of the outflow. The kicked gas particles shock with the gas in the disc immediately after being launched. The net effect is that (i) they slow down, and (ii) they are scattered off their initial trajectory. These effects are more prominent when the particles are launched in a direction closer to the disc plane. Hence, from the case with the outflow along the z -axis, to the case where the direction is slightly tilted (along_jBH0_45 and along_jBH0_10), to the case where the BH spin crosses the disc plane, the maximum outflow speed decreases.
with the gas in the disc immediately after being launched. The net effect is that (i) they slow down, and (ii) they are scattered off their initial trajectory. This effect is more prominent as the particles are launched in a direction closer to the disc plane.
The minor differences in the structure and velocity of the outflows of the three cases presented in this section are consistent with the very small discrepancies observed in the mass accretion rates computed in the three runs, as shown in Figure 6. Overall, some qualitative differences can be seen in the kicked particles distribution up to 300 pc, as illus-trated in the bottom panels of Figure 3. However, when we select the particles according to the radial velocity criterion explained above, those differences become far less evident and the three tests appear similar in their outflow directionθ and their qualitative appearance, up to 400 pc (see Figure 4). The main difference remains instead the outflow maximum speed, that is greater when the particles interact less with the disc. . Mass rates (solid lines) and BH Eddington ratios (dashed lines) for the evolving components of our accretion model, for the four tests described in Sections 3.2 and 3.3 (Ṁ acc−BH computed as in Equation 13,Ṁ in as in Equation 1,Ṁ in after the limit at the self-gravitating mass has been applied, andṀout as in Equation 7). The inflow rate from large scales of run full_model starts to decrease at earlier times, compared to the other tests. This occurs when the BH spin direction crosses the disc plane. Moreover, the sudden increase inṀout (at ∼7 Myr) is caused by the sudden change in the value of the radiative efficiency due to the transition between prograde and retrograde accretion in the sub-grid model (which evolves even when the direction is fixed, as explained in Section 3.1).

The full model
In this section, we present the benchmark run, in which we coupled our sub-grid feedback model with the self-consistent spin evolution model presented in Cenci et al. (2021), as done already for the tests of Section 3.2, additionally modelling an outflow cone axis parallel to the BH spin direction vector, which is evolving with time (as opposed to the cases of the previous section). The initial accretion disc values for this simulation are the same listed in the previous section. The results are shown in Figures 3-6, for an easier comparison with the previous tests.
At around 4 Myr, the BH spin direction starts to be along the CND plane. Since most of the particles are kicked in that direction, the gas is kept far from the BH and, looking at the bottom-right panel of Figure 6, the inflow rate from large scales starts to become lower at earlier times, compared to the other tests. Note that the prescription for the accretion disc evolution limits the disc growth when the inflow rate satisfies our criterion, hence the blue line starts to overlap with the purple one only when the inflow rate is low enough and the disc mass has decreased below the self-gravitating mass.
Furthermore, the large-scale outflow shape is heavily affected by the interaction with the CND, as visible in the bottom-right panel of Figure 3. The polar plot in the last row of Figure 4 (blue line) shows that the outflow is slightly closer to π/4 than in the other three runs. None the less, the outflow direction near the origin, which should be aligned with the launch direction, is already far from it and closer to the z -axis, due to the strong interaction with the disc. Looking at the bottom-right panel of Figure 5, the maximum velocity's magnitude of the outflowing material in this case lowers to ∼300 km s −1 . In this case indeed, the particles, following the BH spin direction, start to be pushed against the disc plane after ∼4 Myr, slowing down even more than in the previous cases.

DISCUSSION AND CONCLUSIONS
We designed a novel BH accretion and momentum-driven feedback model and interfaced it with the publicly available code gizmo. Our prescription is based on that described in Anglés-Alcázar et al. (2017), which we significantly modified. We take into account the presence of a sub-grid accretion disc and ensure self-consistency between its mass evolution, the inflow mass rate from large resolved scales, the outflow mass rate, and the accretion rate on to the central BH. A set of particles is stochastically selected within the BH kernel length, and a fraction of the particles' mass is added to the BH mass. These particles are then kicked (with a velocity magnitude defined by the user) in such a way that the direction is randomly extracted in a cone of a given aperture, whose axis can be chosen arbitrarily, according to the underlying physical assumptions of the model employed. Finally, we coupled the conical feedback model with the accretion disc prescription described in Cenci et al. (2021), wherein the accretion rate on to the BH depends on the inflow rate but is mediated by the accretion disc. Furthermore, we set the reference direction for the conical feedback parallel to the BH spin direction.
There are some caveats in our model that we will address in the future, in the aim at further improving the model. Specifically, (i) when the mass of the sub-grid disc+BH system decreases (due to the decrease in mass of the accretion disc), the BH dynamical mass remains constant, (ii) while statistical properties would be the same, carrying out a new set of simulations and changing the random number generator seeds would produce different positions in the overdensities when the number of particles which exert feedback is small. In addition, particularly in misaligned conditions when the outflow impinges into the CND, the asymmetries could become artificially more pronounced for lower mass resolution in the outflow. Using only the particles which are already present in the simulation to exert feedback restricts the possibilities of implementation to mitigate the above-mentioned effect. A solution could be to generate "wind particles" very close to the BH (to minimise discreteness), which are responsible to create the feedback effect and to which redistribute the dynamical mass in excess (see e.g. Cuadra et al. 2006;Pillepich et al. 2018;Torrey et al. 2020) We carried out a suite of hydrodynamic tests on a CND, comparing the conical feedback prescription with the collimated and radial feedback geometries described in Anglés-Alcázar et al. (2017). We started by exploring how the choice of the cone's axis direction and semi-aperture affects the resulting outflow. Then we performed a test with our complete implementation, using an outflow cone axis parallel to the BH spin direction vector, which evolved with time.
We summarize the findings of our N -body, hydrodynamic tests as follows. Comparing different feedback geometries, we verified that all our tests exhibit the formation of gas overdensities, which interact with the BH and make the accretion behaviour complex (see Figure 2). In the radial case, the particles are heavily scattered when kicked in the direction of the CND plane, and the exerted feedback lowers the accretion considerably with respect to the collimated and conical cases (see Figure 1). In the collimated and conical cases, the kicked particles define and maintain a conical shape.
Analysing the influence of the cone's axis direction and semi-aperture, as well as the full coupling with the BH spin evolution model, we find that: • The shape and orientation of the outflows are similar, regardless of the different orientations and semi-apertures considered in this work (see Figures 3 and 4). As a consequence, also the mass inflows and outflows are similar (see Figure 6).
• A conical shape can always be identified, even at late times, but the outflow's structure is complex (see Figures  3-5). In particular, the outflowing material in the central parts of the outflow has a high speed, which decreases with the distance from the cone's axis (see Figure 5). Moreover, it is not straightforward to relate the launching direction and the resulting outflow orientation (the latter depending also on the distance from the BH), because the interaction with the CND plays a dominant role in shaping the outflow (see Figure 4).
• The more the cone's axis is tilted in the direction of the CND plane, the more the kicked particles are deviated and slowed down (see Figure 5), and the large-scale inflow is decreased (see Figure 6). This is particularly evident when the outflow crosses the CND.
In conclusion, our model is able to self-consistently capture the basic physics of accretion discs and the wind driving mechanism, also accounting for the spin-disc coupling and its effect onto the outflows. It is flexible enough to allow the user to choose a preferential direction for the wind, if desired. This model can be adopted as a sub-grid prescription to multiscale simulations and in future works we plan to employ it to investigate how the spin-feedback coupling affects both the BH growth in mass, the spin evolution, and the interplay between the BH and the ISM of the galaxy host.