{'fraction_non_alphanumeric': 0.03105306749570875, 'fraction_numerical': 0.015268458389849837, 'mean_word_length': 4.525268603263032, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Comparative Study of AGN Feedback Algorithms

Modelling AGN feedback in numerical simulations is both technically and theoretically challenging, with numerous approaches having been published in the literature. We present a study of five distinct approaches to modelling AGN feedback within gravitohydrodynamic simulations of major mergers of Milky Way-sized galaxies. To constrain differences to only be between AGN feedback models, all simulations start from the same initial conditions and use the same star formation algorithm. Most AGN feedback algorithms have five key aspects: black hole accretion rate, energy feedback rate and method, particle accretion algorithm, black hole advection algorithm, and black hole merger algorithm. All models follow different accretion histories, with accretion rates that differ by up to three orders of magnitude at any given time. We consider models with either thermal or kinetic feedback, with the associated energy deposited locally around the black hole. Each feedback algorithm modifies the gas properties near the black hole to different extents. The particle accretion algorithms usually (but not always) maintain good agreement between the mass accreted by \dot{M} dt and the mass of gas particles removed from the simulation. The black hole advection algorithms dampen inappropriate dragging of the black holes by two-body interactions. Advecting the black hole a limited distance based upon local mass distributions has many desirably properties. The black holes merge when given criteria are met, and we find a range of merger times for different criteria. Using the M_{BH}-\sigma relation as a diagnostic of the remnants yields three models that lie within the one-sigma scatter of the observed relation and two that fall below it. The wide variation in accretion behaviours of the models reinforces the fact that there remains much to be learnt about the evolution of galactic nuclei. (abridged)


INTRODUCTION
In the hierarchical model of galaxy formation, the largest galaxies are formed last. Naively, we would expect the highest star formation rates (SFRs) and the activity from Active Galactic Nuclei (AGNs) to occur in these most massive galaxies. However, observational evidence contradicts this, showing that in massive galaxies, the peak SFRs and peak AGN activity occurred at redshifts 1-2 (e.g. Shaver et al. 1996; Madau et al. 1996), and not today. This reduction in activity from z ∼ 2 to today was termed 'downsizing' by Silk & Rees 1998;Ferrarese & Merritt 2000;Gebhardt et al. 2000;Tremaine et al. 2002;King 2003;Gültekin et al. 2009), and the black hole mass and the bulge mass (MBH-M b ; Magorrian et al. 1998;McLure & Dunlop 2002;Marconi & Hunt 2003). A likely implication of these relations is that the AGN feedback is self-regulated: Outflows from the black hole following a strong accretion event interact with the surrounding gas, inhibiting further accretion events, and hence limiting black hole growth (e.g. Silk & Rees 1998;Fabian 1999b;Scannapieco & Oh 2004). Large X-ray cavities in the hot gas halo around AGN (e.g. Boehringer et al. 1993;Mc-Namara et al. 2000) suggest that some outflow events can be quite powerful; these cavities are likely formed by jets from the AGN (e.g. McNamara & Nulsen 2007). Although there is significant evidence supporting the AGN scenario for explaining downsizing, it is not yet understood precisely how feedback energy couples to the surrounding gas; possible mechanisms are line radiation pressure (e.g. Castor et al. 1975;Proga et al. 2000), radiation pressure on dust grains (e.g. Murray et al. 2005), Compton heating of infalling gas (e.g. Ciotti & Ostriker 2001) and photo-ionisation pressure (Buff & McCray 1974;Cowie et al. 1978).
AGN feedback has been implemented in many numerical simulations (e.g. Kawata & Gibson 2005;Di Matteo et al. 2005;Springel et al. 2005;Sijacki et al. 2007;Okamoto et al. 2008;Kurosawa et al. 2009;Booth & Schaye 2009;Debuhr et al. 2011;Dubois et al. 2012), but both theoretical understanding and numerical implementations are fraught with difficulties. First, the spatial scales involved in studying AGN and the related feedback span many orders of magnitude, and information on all of these scales is needed simultaneously to fully and properly understand AGN and their feedback. On the smallest scale is the Schwarzschild radius, where MBH is the mass of the black hole, G is Newton's gravitational constant and c is the speed of light; typical values range rS ≈ 3 × 10 6−8 km. On the next larger scale is the Bondi radius (Bondi 1952), where c∞ is the sound speed of gas at infinity. The value of the Bondi radius is dependent on the black hole's environment and can range from a few parsecs to tens of parsecs (e.g. Kurosawa et al. 2009;Springel et al. 2005). This radius, also known as the capture radius, divides a gas flow in to two regimes (Frank et al. 2002). Consider a spherically symmetric gas cloud centred on a black hole, where the gas is initially at rest at infinity. The only forces acting on the gas are the gravitational force from the central black hole and the pressure forces within the gas (assuming we neglect the self-gravity of the gas). Beyond the Bondi radius, the gas is comparatively uninfluenced by the black hole and flows subsonically. Within r Bondi , the gas density begins to increase, and the gas flow inward eventually reaches a sonic point. At the sonic point, the gas plunges at a free-fall rate in to the black hole (Hobbs et al. 2012). Thus, the Bondi radius can be viewed as the black hole's gravitational radius of influence.
The last spatial scale of interest is that of an entire galaxy (or a galaxy cluster), which can span dozens of kiloparsecs (or a few megaparsecs). When considering all of these scales, comparing the size of a black hole to its massive host galaxy is similar to comparing a coin to the Earth (Fabian 2012).
To complement the range of spatial scales, studying AGN feedback also requires a large range of temporal scales. At short intervals, observations show the luminosity of the central engines of AGN varies on time-scales ranging from days to years (e.g. Webb & Malkan 2000 and references therein; Sarajedini et al. 2011); moreover, there are short term differences in variability amongst the different classes of AGN, making the variable luminosity challenging to understand and constrain. Large scale observations have detected large (∼10 kpc radius) X-ray cavities in the gas around the AGN. To inflate these cavities, an outburst of 10 58 -10 61 ergs of energy would be required (Bîrzan et al. 2004;McNamara et al. 2005) every few ∼10 8 yr, or a timeaveraged output of ∼10 43 -10 45 ergs s −1 (Churazov et al. 2002). Thus, the next inherent difficulty in modelling AGN becomes obvious: AGN luminosity varies on time-scales as short as days, yet they are expected to produce major outbursts every few ∼10 8 yr.
Numerically, we can draw a few parallels between AGN feedback and stellar feedback. Both have been added to numerical simulations to improve realism, and the result was better agreement with observations. However, star formation and stellar feedback is a conceptually simpler problem than AGN feedback. In numerical simulations, if a packet of gas meets a given set of criteria, then stars are formed and feedback energy is returned; the star formation and feedback parameters can be reasonably well constrained using current observations, where there are detailed observations of outflows at many different stages of stellar evolution, such as T-Tauri stars (e.g. Cabrit et al. 1990;Hartigan et al. 1995), Wolf-Rayet stars (e.g. Crowther 2007), mass loss and stellar winds from evolved stars (e.g. Willson 2000) and supernovae (e.g. Reynolds 2008). AGN feedback is much harder to constrain since there are fewer candidates for detailed observation due to the peak of AGN activity being so far in the past. The nearest black hole candidate to study, Sagittarius A*, is currently in a quiescent phase; even during a recent flare to 4 × 10 39 erg s −1 (Nobukawa et al. 2011), its luminosity remained below the typical AGN luminosity range of 10 40 -10 47 erg s −1 (Fabian 1999a). In addition to uncertainty in modelling AGN feedback, tremendous amounts of energy can be numerically returned to the region around the back hole, creating steep gradients and presenting a challenge to numerical integration. Unlike numerical stellar feedback where there is difficulty in preventing unrealistically rapid cooling after the feedback event (e.g. Gerritsen & Icke 1997;Mori et al. 1997;Thacker & Couchman 2000;Springel & Hernquist 2003), continual AGN feedback may create shocks or unrealistically disrupt the system.
When AGN feedback is included in numerical simulations, the approaches can vary widely, in both physics and numerical implementation. Also, the physics is often modelled using different numerical codes with different star formation algorithms and different initial conditions. Even starting from the same initial conditions, mass and spatial resolution, the different hydrodynamical approaches, gas cooling, star formation and stellar feedback algorithms of different numerical codes produce different results, as shown in the comparison by Scannapieco et al. (2012); this comparison did not include AGN feedback. Currently, a direct numerical comparison of the AGN feedback mechanisms does not appear to have been published nor can it be compiled from the assembled literature. Here, we run four of the algorithms found in the literature using the numerical code Hydra. We also present a fifth simulation that is new to this study; this model is designed to take advantage of the 'best' features of the other models. All of our simulations start from the same initial conditions and use the same star formation algorithm; this approach is conservative, but well constrained. We explicitly state that the goal of this is to highlight the different behaviours of the different algorithms and not to critique the various approaches.
By using a pseudo-multiphase star formation model for all simulations, we unavoidably introduce some compromises in the AGN feedback models that rely upon a multiphase gas description; to compensate for this, we have introduced additional variables to represent hot and cold fractions of the gas. Although this star formation model is not precisely equivalent to those implemented in other simulations of AGN feedback, we consider using it a necessary compromise to ensure that the variations in our results are only from the AGN feedback algorithm.
The layout of this paper is as follows. In Section 2 we discuss our simulations, focusing on how the initial conditions are constructed and the AGN feedback algorithms. In Section 3, we will discuss our results, focusing on the impact of the different AGN feedback algorithms. In Section 4, we discuss the final state of each simulation, and we conclude with a review in Section 5.

NUMERICAL SIMULATIONS
To perform our simulations, we use the parallel version of Hydra (Couchman et al. 1995;, which uses an Adaptive Particle-Particle, Particle-Mesh algorithm (Couchman 1991) to calculate gravitational forces and the standard smooth particle hydrodynamics method (SPH; Gingold & Monaghan 1977;Lucy 1977) to calculate gas forces. It includes a star formation algorithm (see section 2.3), and has been modified to include black holes and AGN feedback (see section 2.4).

Galaxy models
To construct our model galaxy, we first use the GalactICs package (Kuijken & Dubinski 1995;Widrow & Dubinski 2005;Widrow et al. 2008) to create a Milky Way-sized galaxy that consists of a stellar bulge, stellar disc, and a dark matter halo; this is done through an iterative process to produce a self consistent system. The free parameters are chosen such that the component masses are similar to the component masses in Springel et al. (2005), and are given in Table 1. The stellar bulge density profile is given bỹ which yields a Sérsic law for the projected density profile if p = 1 − 0.6097/n + 0.05563/n 2 , where n is a free parameter. The constant ρ b is defined using , where σ 2 b is the depth of the gravitational potential associated with the bulge and Re is the radial scale parameter, which is a free parameter; the variable b is adjusted such that Re encloses half the total projected light or mass.
The stellar disc has a truncated density profile that falls off approximately exponentially in R and follows sech 2 in z; the disc has radial and vertical scale heights R d and z d and truncation distances Rtrunc and ztrunc. The radial velocity dispersion is given by σ 2 R (R) = σ 2 R0 e −R/Rσ , where σR 0 is the central velocity dispersion and Rσ = R d for simplicity.
The dark matter halo profile follows where a h is the halo scale length, r h is the cutoff radius, γ is the central cusp strength and σ h is a (line of sight) velocity scale that sets the mass of the halo. The truncation function, C(r; r h , δr h ) = 1 2 erfc r−r h √ 2δr h , smoothly goes from one to zero at r = r h over width δr h .
We then modify this initial galaxy in three ways. First, we convert ten per cent of the total stellar mass into gas to create a gas disc; the gas disc is the same as the stellar disc except that it has been reflected in the x = y plane to avoid coincidence with the star particles. Although the gas scale height is initially larger than physically motivated, cooling allows the gas to collapse into a thin disc within a few 10 Myr; this produces a short transient evolution of the gas accompanied by a brief increase in the SFR. At resolutions higher than presented here, this vertical collapse produces a strong ring-shaped shock which propagates outwards. One solution is to relax the gas disc in a fixed potential and then implement the relaxed disc in the initial conditions; another solution is to reduce the initial scale height of the gas disc. Studies in Williamson & Thacker (2012) show that at our fiducial resolution the pressure wave dissipates comparatively quickly, and is largely insignificant at our lowest resolutions. Given also that we are interested in the relative differences between the feedback algorithms during the merger, which is clearly a much more dynamic event, we consider the above gas disc construction acceptable.
Second, we add a hot gas halo (hgh), which is chosen to follow the observationally motivated β-profile (e.g. Cavaliere & Fusco-Femiano 1976): where ρ0 is the central density, rc is the core radius, and β is the outer slope parameter; we choose rc = 1.75 kpc and β = 2/3 as done in Moster et al. (2011). We set ρ0 by choosing the mass of the hot gas halo within 40 kpc to be equal to two per cent of the total disc mass (Rasmussen et al. 2009); to conserve total halo mass generated by GalactICs, we reduce the mass of the dark matter halo by the total mass of the hot gas halo. By assuming isotropy and hydrostatic equilibrium, the temperature profile of the hot gas halo is given by (Kaufmann et al. 2007) where µ is the mean molecular weight, mp is the proton mass, kB is the Boltzmann constants, and M (r) is the total mass interior to r. The hot gas halo is given an initial angular momentum which scales with the circular velocity, j hgh (R) ∝ Rvcirc(R), where R is the distance from the spin axis of the galaxy (Moster et al. 2011). The final modification to the galaxy is the addition of a black hole (sink) particle, which is placed at the centre of mass and given an internal mass of 10 5 M . While this initial mass is more than ten times lower than expected from estimates of the mass of the Milky Ways central black hole (MBH ∼ (4.36 ± 0.42) × 10 6 M ; Gillessen et al. 2009), it does match the initial black hole masses in Springel et al. (2005) and Debuhr et al. (2011). Moreover, an initial mass lower than anticipated from the MBH-σ relationship would be expected to grow quite quickly in a few Salpeter times because feedback at low masses is comparatively weak. At every step, we calculate a smoothing length for the black hole in order to calculate the gas properties around it, but the sink particle itself only experiences gravitational forces.
Each galaxy has a total mass of 9.60 × 10 11 M , 1 287 743 (168 351) particles for the fiducial (low) resolution simulations, and a Plummer softening length of Plummer = 120 pc ( Plummer = 300 pc). The Plummer softening length is related to the S2 gravitational softening length, S2, by Plummer = S2/2.34. See Table 2 for a breakdown of each galaxy. Lastly, we create a second, identical galaxy, and separate them by 70 kpc; we then place them on parabolic trajectories around one another, in close agreement with the trajectories from Springel et al. (2005).

Verification against other codes
Since many of the models in the literature are run using Gadget2 (Springel et al. 2001;Springel & Hernquist 2002), we ran the same adiabatic merger simulation with both Hydra and the publicly available version of Gadget2 to test the differences. Neither simulation included stellar or AGN feedback, thus we were essentially comparing the gravity and SPH solvers. The galaxies in each simulation followed the same trajectory, and synchronisation -as measured by the time to reach second periapsis -was within 0.2 per cent. Thus we can be confident that our results will be comparable with those currently found in the literature.

Star formation
The star formation algorithm implemented in Hydra is described and tested in Thacker & Couchman (2000); we will provide a brief summary here. This algorithm follows an approach to star formation that has now been studied extensively in the literature, where star formation is allowed to proceed in regions where: (i) The gas exceeds the density limit of nH ∼ 0.01 cm −3 , (ii) The flow is convergent, ∇ · v < 0, (iii) The gas temperature is less than 3 × 10 4 K, (iv) The gas is partially self-gravitating: ρg > 0.4ρDM.
When a gas particle meets the above criterion, the Schmidt Law (e.g. Katz 1992;Kennicutt 1998) is used to determine the amount of gas that is 'converted' in to stars. When the cumulative converted mass of a gas particle reaches half of the particle's original mass, mg, a star particle is spawned with mass mg/2 and the gas particle's mass is reduced to mg/2; when 80 per cent of the remaining gas mass is converted into stars, a second star particle is spawned with mass mg/2 and the gas particle is removed from the simulation. For computational efficiency, Hydra uses a Lagrangian form of the Schmidt Law, based upon the local density of the SPH particle, where C sfr is the star formation rate normalisation, ρg is the SPH density of the gas particle, and Mg is the mass of the gas in the particle that has not been converted in to stars. All other calculations in the code (e.g. density, smoothing length, etc...) assume that the entire particle is a gas particle of mass mg or mg/2. Whenever a star particle is created, feedback energy is immediately returned to the surrounding environment. Although our fiducial resolution models takes O(10 2 ) steps to evolve through the lifetime of an 8 M star, the lack of delay between star particle formation and feedback energy release is not a significant issue because at the resolution considered here we are still averaging over a number of giant molecular clouds per particle and the evolution is representative at best. For every 100 M of stars formed, there is one supernova event, which contributes 10 51 ergs to the interstellar medium (Sommer-Larsen et al. 1999), and feeds back 5 × 10 15 e * ergs g −1 of gas converted in stars, where e * is a dimensionless parameter; we set e * = 0.4 to match Navarro & White (1993).  As first showed by Katz (1992), feedback energy returned to the interstellar medium (ISM) is radiated away quickly in high density regions. A reduced density is used in the cooling algorithm to prevent this immediate loss, thus allowing the feedback to influence the surrounding environment; the reduced density decays back to its local SPH value with a half-life of t 1/2 = 5 Myr. The parameters e * and t 1/2 are chosen to match those set in Thacker & Couchman (2000) to reproduce the Milky Way's star formation rate in a simulation of an isolated Milky Way-like galaxy.
This star formation algorithm is conceptually similar to the methods used in Brook et al. (2004), Stinson et al. (2006) and Governato et al. (2007), however, it varies from the subresolution multiphase models found in Springel & Hernquist (2003), Springel et al. (2005), Booth &Schaye (2009) andHopkins et al. (2009). In multiphase models, the subgrid physics is derived from a model of cold clouds (where stars form) embedded in a hot, pressure-confining phase. Above a given density threshold, the gas is thermally unstable to the onset of this two-phase medium. The mass fraction in each phase is determined by star formation and feedback, evaporation of the cold clouds through thermal conduction, and the growth of clouds through radiative cooling. The star formation rate is then calculated based upon a prescribed law and an equation of state. A noteworthy difference between this method and our method is that here, gas can freely flow between phases, whereas in Hydra, once gas is 'converted' to stars, it is carried forward and not allowed to cool until the specified cooling period is reached.

Black hole and AGN feedback algorithms
AGN feedback algorithms essentially have five key components: (i) The accretion rate on to the black hole, (ii) The SPH particle accretion algorithm, (iii) The energy feedback algorithm, (iv) The black hole advection algorithm, and (v) The black hole merger algorithm.
Each component will briefly be discussed below, then we will discuss the models in section 2.5.

Accretion rates
A commonly used accretion rate is the Bondi accretion rate (Bondi 1952) where ρ∞ and c∞ are the gas density and sound speed at infinity, v is the relative velocity between the gas at infinity and the black hole, and MBH is the mass of the black hole.
In hydrostatic equilibrium, the maximum physical accretion rate is the Eddington accretion rate, where mp is the proton mass, σT is the Thomson cross section, and r is the radiative efficiency (i.e. the mass-to-energy conversion efficiency); we set r = 0.1 (Shakura & Sunyaev 1973). Typically, numerical simulations limit their accretion rate to the Eddington accretion rate, and we follow this convention. If the assumption of spherical symmetry is ignored, super-Eddington accretion rates can be achieved, which are associated with collimated outflows. In small scale outflow simulations by Kurosawa et al. (2009), they find steady-state results where super-Eddington accretion occurs in the equatorial plane.

Black hole mass growth and particle accretion
In all cases, the 'internal' and 'dynamical' masses of the black hole is tracked. The internal mass, MBH, is the mass of the black hole, which is increased byṀBHdt at every iteration; this mass is used in all calculations concerning AGN feedback. The dynamical mass, mBH, is the mass of the sink particle, which is increased by the mass of a gas particle whenever one is accreted. The particle accretion algorithm should ideally maintain MBH ∼ mBH and directly address the loss of gas near the black hole due to accretion. In Model DQM (see §2.5.4), the dynamical mass is fixed for all time (except during a black hole merger), so the included particle accretion algorithm only simulates the loss of gas. When a particle is accreted on to the black hole, its mass and momentum are added to the black hole particle (except for Model DQM), and the gas particle is removed from all further computations; this accretion does not affect the internal mass of the black hole. There are three categories of particle accretion algorithms: (i) Stochastic-Unconditional : At all times, nearby particles are tested to see if they will stochastically accrete.
(ii) Stochastic-Conditional : Nearby particles are tested to see if they stochastically accrete if given criteria of the mBH-MBH relationship are satisfied.
(iii) Continual-Conditional : Nearby particles are continually accreted while given criteria of the mBH-MBH relationship are satisfied.

Feedback
In all simulation, a portion of the accreted mass is returned as feedback energy, E = Ṁ BHc 2 dt, where is a dimensionless efficiency parameter. This energy is returned to the nearby particles (except Model ONB; see §2.5.3) either by increasing their internal energy (thermal feedback models), or by increasing their momentum via p = E/c (kinetic feedback models).

Black hole advection
Properly tracking the black hole particle is critical since accretion rates (hence feedback rates) depend on the local gas properties around the black hole. When the black hole mass is similar to a star or gas particle mass, it can be easily and inappropriately dragged around by two-body forces leading to the inaccurate calculation of gas properties. Avoiding such behaviour is clearly desirable. In all of the simulations, a black hole advection algorithm is implemented to minimise the inappropriate movement; except in Model DQM (see §2.5.4), the black hole's position is artificially updated after the completion of the gravitational solver algorithm but before the calculation of the accretion rate.

Black hole mergers
The process by which black holes merge is still a matter of active research (e.g. Escala et al. 2004;Berentzen et al. 2009;Khan et al. 2011;Bode et al. 2012). Given that our simulations result in a combined halo and stellar system producing drag, it is reasonable to expect that the black holes would merge. Since any merger occurs on sub-resolution scales, models include a merger prescription to instantly merge the black holes when given criteria are met. When two black holes merge, one sink particle is removed and the remaining sink particle has the combined mass (both internal and dynamical) and momentum, and is repositioned to the centre of mass of the two progenitors.

Black hole's local environment
Around every black hole we define a radius of influence, r inf . All gas particles within r inf contribute to the accretion properties at the black hole, are eligible to receive feedback energy (except for Model ONB), and are eligible to be accreted on to the black hole particle. Gas particles outside of r inf have no explicit impact on the black hole or the AGN feedback algorithms. We set r inf = max(2hBH, 2hmin), where a sphere with radius 2hBH around the black hole particle includes 60 gas particles, and hmin is the smallest resolved smoothing length in the SPH solver.

The Models
In sections 2.5.1 to 2.5.5, we describe the particular AGN feedback algorithms of our five primary models; a summary of each model can be found in Table 3. We ran four additional models, each of which is a slight variant of a primary model; these are described in section 2.5.6. Every model was run at both fiducial and low resolutions. For nomenclature, we name these models after the initials of the authors of the paper it originally appeared in. The variant models have the same name as their parent model, followed by a lower case character to signify the difference. A model name followed by a subscript 'l' explicitly refers to the low resolution version of the model; no subscript will refer to either the fiducial resolution version or to both versions, depending on context.

Model SDH
This model is based upon that found in Springel et al. (2005) (herein SDH05). The accretion rate is given by a modified Bondi accretion rate, where cs and ρ are the local sound speed and density of the gas, and v rel is the relative velocity of the black hole to the nearby gas. The free parameter, α, is included to relate the numerically calculated gas density and sound speed to what one would expect in reality. Booth & Schaye (2009) argue that modest resolutions underestimate the density and overestimate the sound speed by orders of magnitude, thus justifying large values of α. As in SDH05, we set α = 100. Finally, the accretion rate is Eddington-limited, thusṀBH = min Ṁ B,ṀEdd . A given fraction of the accreted mass is allowed to return to the surrounding environment as feedback energy. The rate of return isĖ where r = 0.1 is the radiative efficiency, and f = 0.05 is the fraction of energy that can couple with the gas. This energy is returned kernel-weighted to the gas particles within r inf .
To track the growth of the black hole particles, a stochastic-unconditional particle accretion algorithm is used. At every iteration, a probability, pi, is calculated for each particle, i, within r inf . This probability is then compared to a random number, xi, and if pi > xi, then the particle is accreted. In this algorithm, the probability is given by where wi is the kernel weight of gas particle i relative to the black hole, and xi is uniformly distributed on the the interval (0,1).
To minimize inappropriate motions of the black hole particles, at every iteration the black hole is relocated to the nearby gas particle with the minimum potential energy provided that v rel < 0.25cs. If no gas particle meeting the velocity criteria exists within r inf , then the black hole is not advected. Once mBH > 10mg, the black hole advection is Particle Accretion stochastic-unconditional stochastic-conditional continual-conditional continual-conditional continual-conditional BH Advection gas particle with gas particle with ∆l ONB tracer mass ∆l WT v rel < 0.25cs with v rel < 0.25cs along stellar gradients towards centre of mass and smallest U and smallest U Table 3. A summary of our five primary models describing each key algorithm. More detail and definitions of the variables are given in subsections 2.5.1 to 2.5.5. Numerically, h BH = max(h BH , h min ), where a sphere with radius 2h BH around the black hole particle includes 60 gas particles, and h min is the smallest resolved smoothing length in the SPH solver.

BH Merger
turned off and its movement is handled only by the gravitational solver. Lastly, two black holes merge when they come within each other's smoothing lengths and have a relative velocity less than the local sound speed; SDH05 argue that the local sound speed represents a simple measure of the velocity scale at which the black holes are able to merge.

Model BS
This model is based upon that found in Booth & Schaye (2009) (herein BS09); this model has many similarities to Model SDH, but was originally implemented in a cosmological simulation. We caution that implementing their model in a higher-resolution simulation could lead to unwanted behaviours as a number of the model parameters were chosen for their specific resolution. This model uses the (modified) Bondi accretion rate given in (10). If the resolution is sufficient enough, BS09 argue that the justifications for α = 100 given in section 2.5.1 break down in low-density regions. Thus, this model sets α to be a function of the local hydrogen density, nH: where n * H is the critical value required for the formation of a cold interstellar gas phase, and β is a free parameter; as in BS09, we set n * H = 0.1 cm −3 and β = 1. Finally, the accretion rate is Eddington-limited, thusṀBH = min Ṁ B,ṀEdd .
The rate of feedback is also given by (11), with r = 0.1 and f = 0.15, which is three times more efficient than Model SDH. The feedback energy is allowed to accumulate until E > Ecrit, at which point a random gas particle within r inf receives Ecrit energy; this is repeated until the accumulated energy drops below Ecrit. The critical energy is defined as where mg is the (initial) mass of a gas particle and ∆T is the temperature increase a particle experiences with every feedback event. Due to the higher resolution in our model (m gas, fiducial = 3.6 × 10 4 M compared to 8.64 × 10 7 M h −1 in BS09; fiducial = 120 pc compared to 2 kpc h −1 after z = 2.91 in BS09), we set ∆T = 5 × 10 6 K as opposed to ∆T = 10 8 K used in BS09. Using the value in BS09 produces very large temperature gradients in the SPH solver that are difficult to integrate accurately. Gas particles are accreted by a stochastic-conditional particle accretion algorithm. If MBH < mBH, then the probability of accretion is pi ≡ 0, otherwise it is calculated using As in Model SDH, particle i is accreted if pi > xi.
The black hole advection is the same as in Model SDH. Lastly, two black holes merge when they come within each other's smoothing lengths and have a relative velocity less than the circular velocity at the radius of the most massive black hole's smoothing length. BS09 explicitly state that this merging criteria differs from SDH05 since the feedback returned may temporarily increase the local sound speed, thus may not be a representative velocity scale.

Model ONB
This model is based upon the model found in Okamoto et al. (2008) (herein ONB08), which was originally implemented in a simulation with cosmological initial conditions. We also note that this model is distinctly different from our other models in that it is specifically designed to reproduce the radio mode of feedback. In this model, it is assumed that radiation from stars (e.g. through starbursts or winds) interacts with the rotating, clumpy ISM. This radiation irradiates one layer of gas at a time, extracting angular momentum from it. This permits an inflow of gas towards the galactic centre, and ultimately on to the black hole itself (Umemura et al. 1997;Umemura 2001;Kawakatu & Umemura 2002). Thus, an accretion rate of gas on to the black hole can be calculated by considering the stellar clouds in the region of star formation (RSF) near the black hole. Using these assumptions, ONB08 calculate the drag accretion to bė where drag = 1 is the drag efficiency, LRSF is the total bolometric luminosity of all the stars in the RSF, and τRSF is the total optical depth of the RSF. The total luminosity is calculated by summing the age-dependent bolometric luminosities, which are obtained from a lookup table generated by pegase2 (Fioc & Rocca-Volmerange 1997). Next, the total optical depth of a gas cloud is given by τc = χ d ρcrc χ d mc/r 2 c where ρc, mc and rc are the density, mass and radius of the cloud, respectively. Assuming that all the clouds are identical and randomly distributed over the region of star formation, the total optical depth can be approximated by where χ d = 50 cm 2 g −1 is the mass extinction coefficient, Mc is the total mass of the clouds in the RSF, and RRSF is the radius of the RSF; since Hydra does not explicitly track a multiphase gas, we set Mc to be half of the total gas mass within the RSF. We initialise RRSF = max (R40, 2hmin), where a sphere with radius R40 centred on the black hole contains 40 gas particles. Then, larger radii are searched to find the radius that maximisesṀ drag . Lastly, if the gas density within the sphere is less than ρ thresh = 5×10 −25 g cm −3 , the accretion rate is set to zero.
In this model, it is explicitly assumed that the feedback heats the halo gas through the production of jets. The jet mechanism used here is that of Meier (2001), and generates power from the rotational energy of the accretion flow and from the black hole itself. The accretion flow is divided into two regimes: standard thin discs (SD; optically thick, geometrically thin, radiatively efficient) and radiatively inefficient accretion flows (RAIF; optically thin, geometrically thick). Using the parameters in ONB08, the respective accretion-dependent luminosities are whereṁ ≡ṀBH/Ṁ Edd andṁcrit ≈ α 2 is the critical accretion rate that sets the division between the SD and RAIF regimes (e.g. Narayan et al. 1998). As in ONB08, we set αSD = αRAIF ≡ α = 0.1, where αSD and αRAIF are the viscosity parameters for SD and RAIF's, respectively. The feedback rate is then given byĖ feed = rLjet, and the energy is distributed equally to the 40 nearest diffuse gas particles with ρ < 0.1ρ thresh .
To accrete particles, this model uses a continualconditional particle accretion algorithm; when the internal black hole mass exceeds its dynamical mass, nearby gas particles are accreted with the probability where NRSF is the number of gas particles within the region of star formation, and mi is the mass of the i'th gas particle; particles are accreted until the dynamical mass exceeds the internal mass.
To track the black hole, at every iteration, the local stellar density fields are computed and the black hole is moved along the steepest gradient by an amount ∆lONB = min(0.01 S2, 0.03 |v| dt), where S2 is the gravitational softening length, |v| is the velocity of the black hole, and dt is the time-step; these coefficients are the same as in ONB08 and were determined empirically.
Lastly, two black holes are assumed to merge when they are within a softening length of one another and are gravitationally bound.

Model DQM
This model is based upon that found in Debuhr et al. (2011) (herein DQM11), and models accretion through the transportation of angular momentum, which is based upon multiscale SPH simulations by Hopkins & Quataert (2010). The accretion rate isṀ where δ is the dimensionless viscosity, Σ is the mean gas surface density, and Ω = GM/r 3 inf is the rotational angular velocity of the gas. DQM11 treat δ as a free parameter, and it characterises both the efficiency of angular momentum transport and the fraction of gas that is being converted into stars versus being accreted on to the black hole; as in DQM11, we set δ = 0.05.
The feedback is returned as momentum usinġ where τ is the infrared optical depth, and L = min rṀviscc 2 , L Edd is the Eddington-limited luminosity.
This momentum outflow is used to approximate radiation pressure produced by absorption and scattering of the AGN's feedback; specifically, the ultraviolet radiation (emitted from the black hole) will depositṗUV = L/c on to the gas, while infrared radiation (re-emitted from dust) will de-positṗIR = τ L/c on to the gas. Thus, the totalṗ can be approximated as (1 + τ )L/c τ L/c for τ 1, which is valid near the peak of AGN activity when the black hole gains most of its mass. As in DQM11, we set τ = 10. Lastly, the momentum is returned radially, such that every gas particle within r inf receives an equal acceleration. In this model, there is no explicit artificial black hole advection algorithm. Instead, a tracer mass is used to represent the black hole particle; that is, the dynamical mass of the black hole is initialised to mBH = 10 9 M and held constant throughout the simulation. The black hole particle is now only advected by the gravitational solver; since the tracer mass in the fiducial resolution simulations is ∼10 3 mDM and ∼3 × 10 4 mg, it will not undergo artificial dragging by the surrounding particles. By initialising mBH = 10 9 M rather than mBH = MBH, the mass of each galaxy is increased by 0.01 per cent. This means that Model DQM has slightly (but unavoidably) different initial conditions than the remainder of the models. The internal mass of the black hole is still initialised to 10 5 M , and this mass allowed to grow as calculated by the accretion rate given in (22).
Since the dynamical mass of the black hole is fixed, a continual-conditional particle accretion algorithm removes 'accreted' gas particles, but does not add their properties to the dynamical mass. Here, a random particle within r inf is removed whenever there is a mismatch between the amount of gas accreted via (22) and the total mass of removed gas particles.
Lastly, black holes merge when they come within r inf of each other, regardless of velocity. In this model, black hole mergers are the only mechanism to increase the dynamical mass of the black hole.

Model WT: This study
This model uses the modified Bondi accretion rate given in (10) for both its simplicity and for its wide use in the literature (e.g. SDH05; Robertson et al. 2006;Croft et al. 2009;Johansson et al. 2009); as in Model SDH, we set α = 100. The energy feedback rate is given by (11), and the feedback energy is equally distributed amongst all the gas particles within r inf . This provides an isotropic heating to the core rather than preferentially heating (or super-heating) the particles very near the black hole, which may be transient.
The black hole advection algorithm is a modified version of that presented in Model ONB. First, the black hole is displaced towards the centre of mass of the sphere with radius r inf centred on the black hole rather than along stellar gradients. This method still gives preference to the stellar distribution, but also considers the dark matter and gas distributions. Second, the distance the black hole is displaced has been modified to ∆lWT = min(0.10hBH, 0.30 |v| dt, dCM), where dCM is the distance from the black hole to the centre of mass. We choose to use the black hole smoothing length rather than the softening length since all of the properties near the black hole are calculated using 2hBH as the characteristic length. As in (21), our coefficients are empirically chosen so that non-negligible displacement would be possible. Specifically, for the first option, hBH < S2, so its coefficient needs to be increased; for the second option, our resolution is higher than in ONB08, thus we will have a smaller dt, thus we again need a larger coefficient. We choose this method rather than potential-well method of Models SDH and BS for the following reasons: If the AGN feedback creates a gas void around the black hole, then either there will be no artificial displacement (allowing for the possibility of unnatural motion) or the black hole will be coupled to a gas particle on the edge of the void. In this method, as in Model ONB, the black hole only moves towards a particle rather than being coupled to it. As in Models ONB and DQM, we include a continualconditional particle accretion algorithm: When MBH > mBH + mg/2, we accrete the gas particle that is nearest to the black hole. The term mg/2 forces the dynamical and internal mass to oscillate around one another, and choosing the nearest particle is to reduce random effects in the simulation.
The black hole merger algorithm is the same as given in Model SDH. We choose this algorithm since it had a distance and a velocity criteria, and its velocity requirement is less stringent than in Models BS and ONB. The merger algorithms in Models BS and ONB were created for cosmological models, whose resolution is lower than our fiducial resolution runs, thus we argue that our implementation of their merger prescriptions are more stringent than they intended.

Additional models
We have tested four additional models, each of which is a slight variation of one the above models. These models are identical to their parent models described above, except for the variation listed below.
• Model BSw : This uses the black hole advection algorithm of Model WT; early tests of Model BS showed very erratic black hole motion that could possibly compromise the results.
• Model ONBc: This uses a very conservative search algorithm to calculateṀ drag ; this yields a very small RRSF and a very low accretion rate.
• Model WTh: For a resolution test, this model uses hmin → hmin/2; this should only impact calculations performed in very dense regions.
• Model DQMe: This model uses r inf ≡ 4 S2 = 1.17 kpc (compared to r inf ∼ 73 pc of Model DQM); in DQM11, they fix r inf ≡ 4 S2, although this value is ∼188 pc in their models.

RESULTS
Each of our models was evolved through a merger event, similar to that of SDH05. Model ONB was evolved for 1.1 Gyr and the rest were evolved for 1.5 Gyr (including Model ONB l ). By returning the feedback energy to the halo gas in Model ONB, a dense galactic core formed near the core merger epoch, which resulted in an extremely large wallclock time per step for the fiducial resolution version.
Each model followed a similar qualitative history, which is shown for Model WT in Figs. 1 and 2 for the gas column density and gas temperature, respectively. From this evolution, we see four significant epochs: first periapsis, apoapsis, second periapsis, and core merger. The times of these epochs are calculated using the black holes as proxies for the centres of the galaxies, thus the epochs essentially represent the local minimum and maximum separation of the two black holes. In all models, first periapsis occurs at 166 Myr and apoapsis occurs at 480 Myr; the latter is sustained for approximately 100 Myr. Second periapsis occurs approximately 11 Myr earlier for Model DQM than for the rest of the models. We have verified that this more rapid evolution is due to the inclusion of a more massive sink particle representing the black hole. In two-body simulations, where each particle represents one galaxy and dynamical friction is avoided, the period of the particles representing galaxies in Model DQM is 3.4 Myr (0.23 per cent) shorter than the particles representing the galaxies without tracer masses. For the full system, we find that the maximum separation at apoapsis is very slightly smaller (1.1 per cent), likely due to slightly higher dynamical friction, leading to an earlier second periapsis. Lastly, the core merger takes approximately 10 Myr to complete, and the onset of this process spans a range of 15 Myr amongst the different models. See Table 4 for a list of when these epochs occur.
As might be expected, there are notable morphological differences between the models; this is readily apparent in Fig. 3, which displays the gas temperature, gas column   Figure 4 where we have plotted the (S2 softened) rotation curve for five variations on Model WT l just prior to the onset of bar formation; in this test suite, a tracer mass represents the black hole and has been varied between 10 6 and 10 9 M . The rotation curves for Mtracer = 10 6−8 M are similar, and the peak velocity is 2.5 per cent higher for Mtracer = 5 × 10 8 M and 4.0 per cent higher for Mtracer = 10 9 M . The higher rotational velocities for Mtracer = 10 9 M stabalises the galaxy against bar formation.
To ensure that our initial black hole mass is reasonable, we ran Model WT l with seed black hole masses of 10 5 , 10 6 and 10 7 M . All three models had final black hole masses within 1.5 per cent of one another, and all three remnants had MBH-σ relations that agreed with the observed relation within the one-sigma standard deviation.

Black hole advection
The location of the black hole with respect to the galactic core plays a fundamental roll in determining the gas properties used to calculate the accretion and energy feedback rates. Even small displacements with respect to the centre of the potential well can have notable effects. We have studied four artificial advection algorithms, each yielding different results. Since the results presented in the rest of this paper are implicitly coupled to the behaviour of the advection algorithm, we discuss this issue first. See Fig. 5, where we have isolated the path of one black hole in Models SDH, ONB and DQM.
Without any artificial advection, low mass black hole particles are easily dragged by the local particles (whose mass may be comparable or greater than the mass of the black hole); even with artificial advection, some black hole 'chaotic' motions (or oscillations) persist, depending on the algorithm. The smoothest path is from the tracer mass (Model DQM); however, the increased mass of the galaxy slightly adjusts its trajectory, decreases its evolution time compared to the rest of the models, and impacts the disc morphology (i.e. prevents bar formation). The least smooth path is from coupling the black hole to the gas particle with the minimum potential energy, provided that v rel < 0.25cs (Models SDH and BS). Frequently, there is no gas particle  satisfying the velocity criterion, thus the black hole oscillates about the galactic core. When there is a gas particle satisfying the criteria, it is likely near the periphery of r inf since the inner particles are typically too hot. Thus, the stringent velocity criterion permits both dragging by two-body forces and large artificial jumps in position. In Models SDH and BS, the artificial advection is turned off when mBH > 10mg; in the fiducial (low) resolution, the motions become smaller (larger) and more chaotic. We note that SDH05 and BS09 both use a multiphase star formation model, which produces a smoother temperature profile around the black hole than our star formation model. This difference may explain the poor advection in Models SDH and BS. From the simulations, the advection methods in Models ONB and WT appear to have some advantages; they do not add mass to the galaxy as in the tracer mass method, nor do they couple to gas particles. By associating the direction of motion to stellar densities (Model ONB) or the local centre of mass within r inf (Model WT), undo influence of gas particles, which are dynamically and continuously affected by the AGN feedback, is avoided. Lastly, by limiting the distance a black hole can be artificially advected per iteration, it is possible for it to remain in the centre of a void.
The advection algorithms in Models DQM, ONB and WT are comparatively unaffected by resolution; this is expected since the tracer mass is independent of the local gas particles, and local averages vary little with resolution. The advection algorithm in Models SDH and BS are more resolution dependent; with lower resolution comes a larger r inf , thus the black hole has a larger distance it can be artificially displaced. It is worth noting that Fig. 5 clearly shows that the oscillations of Model SDH is damped with higher resolution, but is not completely removed.
Due to the notable black hole motion in Model BS, we ran a variation, Model BSw, which uses the black hole advection algorithm of Model WT. In this model, the black hole had smoother accretion rates and more continuous gas densities and temperatures near the black hole. The outbursts near apoapsis were stronger and more symmetric between the two galaxies than in Model BS. Ironically, although we might argue that Model BSw is inherently better than Model BS because of the better black hole trajectories, both models yield remnants with similar total black hole masses and stellar velocity dispersions, even though the black holes in Model BS never merge and the feedback history of the two models is different.

Accretion rate
The mass evolution of the black hole is a key aspect of the AGN feedback models and is an important diagnostic since it sets the MBH-σ results. In the top panel of Fig. 6, we show the total black hole mass over time; in the bottom panel, we show the accretion rates (geometrically averaged over both black holes and plotted in bins of 10 Myr); in the middle panel, we show when a black hole is accreting at its Eddington limit. There are three major epochs of black hole growth: at the beginning (Models ONB, WT and DQM), at apoapsis (Models SDH, BS and WT), and at core merger (Models SDH, BS and DQM). It is typically at these times when a black hole is undergoing Eddington accretion; see Table 5 for the percent of time that a black hole is accreting at its Eddington limit.
Three models use the Eddington-limited Bondi accretion rate: Models SDH, BS and WT. In SDH05 (see bottom panel of their fig. 14), the accretion rate increases from the onset of the simulation until prior to second periapsis, at which point there is only a slight decrease. At core merger, there is an increase of approximately 1.8 dex, followed by a sudden drop as the system is totally disrupted. Since our initial conditions and our Model SDH was constructed to mimic SDH05, we expect similar results, which, within reason, are obtained. In both Models SDH and WT, there is an initial decline in accretion rate as the galaxies are relaxing; this initial relaxation does not occur in SDH05 since the vertical structure of the gas is set in hydrostatic equi- librium, thus this initial difference is expected. After first periapsis, gas is funnelled into the core, thus, in agreement with SDH05, there is an increasing accretion rate. At apoapsis, there is a plateau in Model SDH, but a short-lived spike in Model WT. At core merger, there is an increase followed by a rapid decrease as the systems are totally disrupted. For the first 400 Myr, the accretion rate of Model BS follows the same trends as Models SDH and WT. During this time, α(nH) > 100, which explains whyṀ Model BS > M Models SDH & WT for 230 t/Myr 400. After this point, M Model BS decreases, which is a result of the strong outflows from feedback and the motion of the black holes; the latter prevents the black holes from remaining in a steady environment, resulting in rapidly varying gas characteristics around the black hole. The accretion rates in Model BSw continue to increase until apoapsis, after which they decrease due to strong outflows, but to a smoother accretion profile that is slightly higher than in Model BS; see Fig. 7.
In DQM11, there is a large jump inṀvisc followed by a slow decline shortly after first periapsis and again starting at second periapsis; see the bottom panel of their fig.  1. Their initial galactic separation is approximately twice of ours, thus their evolution times are longer, which may account for differences between our models and theirs. Also, their mass and spatial resolution are slightly better than in the models presented here, thus their feedback is likely distributed amongst fewer particles, which keeps their void relatively smaller (see discussion in §3.3).
In Model DQM, there is not the accretion epoch between first periapsis and apoapsis that occurs in DQM11; this is a result of both a shorter time during which gas can fall into the galactic cores, and the slightly larger r inf which allows for efficient feedback to maintain the void during non-cataclysmic events. At second periapsis, there is an increase of approximately 2.2 dex in Model DQM compared to DQM11's increase of approximately 3 dex. After the merger, both accretion rates drop off, but DQM11 falls off faster.
ONB08 is a cosmological simulation run from z = 49 to z = 0. Over the final 7 Gyr of their simulation, there are no significant galaxy mergers and there is negligible black hole growth (see their fig. 8). During this time, their accretion rate lies comfortable within the RAIF regime by more than an order of magnitude. In Model ONBc, there is a low accretion rate that is continually decreasing and comfortably within the RAIF regime; however, this model uses a very conservative search algorithm to calculateṀ drag , thus the low accretion rate is likely artificially low because of the numerical scheme. When we use a less stringent search algorithm that is better suited for use in Hydra at the resolutions we present, the accretion rate quickly rises tȯ M drag ∼ 0.005 M yr −1 and remains approximately constant for the remainder of the simulation. This relatively constant rate is a result of the feedback energy being deposited into the halo gas rather than modifying the galactic core. Since the feedback is being returned to the halo gas, there are negligible core differences between Models ONB and ONBc, except for the accretion rate.

Energy feedback algorithm
The two broad categories of returning the feedback energy are thermally,Ė, and kinetically,ṗ. Thermal energy is returned to gas particles in Models SDH, BS, ONB and WT by increasing the internal energy of the particles; kinetic energy is returned to gas particles in Model DQM by increasing the momentum of the particles. Although these two approaches are related byĖ =ṗc, both forms of feedback yield drastically different results, which are readily apparent on both large scales and small scales (e.g. see Fig. 8 for gas density and temperatures around the black holes).
The feedback in Model ONB is returned to the gas particles with ρ < 5 × 10 −26 g cm −3 ; thus, this energy is distributed to halo gas particles with distances 0.5 r/kpc 3.5 (7 r/r inf 48) from the black holes, which leaves the galactic cores relatively unmodified by AGN feedback. As gas falls towards the galactic cores (due to tidal interactions at first periapsis), it remains there since the AGN feedback is not local enough to remove it, hence these models have the highest core density. The major core heating events in this case are from shock heating of the infalling gas after first periapsis and from the final core collisions. The AGN feedback in the remaining models is delivered directly to the core region, leading to a more distinct modification of the entire system. From the onset of the remaining thermal feedback simulations, there is a higher core temperature and a lower core density than in Model ONB. Tidal interactions at first periapsis cause gas to be funnelled towards the cores, but unlike Model ONB, the temperature continues to rise from AGN heating, and the density remains moderated. Shortly after first periapsis, a galactic bar starts to form and persists strongly for approximately 200 Myr; as it dissipates there is a nuclear inflow of gas, accounting for the local density maxima near apoapsis in Models SDH and WT. After the dissipation of the bar, the galaxies begin to eject gas, preferentially along the polar axis. As the cores merge, there is a final set of feedback events as the gas is heated and blown away; after this, accretion rates and core densities drop while core temperature remains approximately constant.
Kinetic feedback affects the galactic core in different quantitative and qualitative ways than thermal feedback. In this model, the energy is returned radially to gas particles as a momentum boost, which quickly and efficiently forms a void around the black holes; in Model DQM, a void of r ∼ 0.55r inf is formed within the first 60 Myr. Once the void is formed, it persists for the remainder of the simulation. Strong events, such as the infalling gas after first periapsis or the collision at core merger, can briefly decrease the void radius, increasing the core density. Otherwise, the void and core characteristics are held constant, as we can see by an approximately constant temperature and density profile between first periapsis and core merger.
Model DQMe uses the same kinetic feedback algorithm as in Model DQM, but instead r inf ≡ 4 S2, which leads to an initially higher accretion rate. The increased amount of feedback energy is returned to more particles, delaying the formation of the void. When the void forms after 140 Myr, it has a radius of r ∼ 0.85r inf ∼ 1 kpc, which is highly unrealistic. This void is more rigid than the void in Model DQM, only being affected at second periapsis and core merger. The unrealistically large r inf in Model DQMe also leads to a premature merger of the black holes, which is poorly handled by the merger subroutine. Thus, the importance of choosing a physical r inf is clear: catastrophic and unphysical results will ensue if it is not chosen appropriately.

Particle accretion algorithm
A particle accretion algorithm should accrete particles such that the internal mass, MBH, and dynamical mass, mBH, of any given black hole remain similar for all time. In Model DQM, a tracer mass is used to represent the black hole particle, rendering the comparison between internal and dynamic masses moot; these models use a continual-conditional algorithm to remove 'accreted' gas, thus any conclusions we make regarding this category can be applied. In the remainder of our fiducial resolution models, the first accretion event will increase the dynamical mass by 36 per cent since mBH = 0.36mg. All subsequent accretion events will add relatively less mass to the black hole particle. Thus we expect the ratio ∆m/MBH, where ∆m = mBH − MBH, to start at |∆m/MBH| 0.36 and decay to zero with time. In the top panel of Fig. 9, we plot ∆m/MBH for one black hole in Models SDH, BS, ONB and WT. The expected behaviour is obtained in Models BS, ONB and WT; the decline in Models ONB and WT is smoother than in Model BS since there is no stochastic component to the accretion. The expected decline is not observed in Model SDH, where the ratio has no apparent trend, and at late times ∆m/MBH → −0.01 rather than zero.
The black holes do not all grow at the same rate. Thus, for a normalised comparison, we plot ∆m/mg in the bottom four panels of Fig. 9 for Models SDH, BS, ONB and WT, respectively. Since every accretion event adds mass mg or mg/2 to the dynamical mass, we expect |∆m/mg| 1 for all time, unless multiple particles are accreted in one step, which can occur in the stochastic models.
In Model SDH, particle accretion events are only dependent on the environment. This lack of MBH-mBH coupling permits the internal and dynamic masses to follow different histories, as seen by the lack of trend in the second panel of Fig. 9. The discontinuity in ∆m/mg at core merger is a result of both black holes following different mass growth histories, both internally and dynamically. The stochastic-conditional accretion algorithm in Model BS allows for fluctuations in an otherwise smooth process, and the probability's dependence on ∆m forces good agreement between the internal and dynamic masses. Finally, by design, Model ONB (WT) always maintains 0.0 ∆m/mg 1.0 (−0.5 ∆m/mg 0.5), displaying the effectiveness of the continual-conditional particle accretion methods.  Table 6. Merger time of the black holes, as measured from the beginning of the simulation. Dashed lines indicate that the black holes did not merged within our 1.5 Gyr simulation time.

Black hole merger algorithm
All of the models include a black hole merger prescription, which conserves internal, dynamical and tracer mass; see Table 6 for precise merger times. Thus, the merger should have limited gravitational impact on the surrounding environment. However, the amount of feedback energy returned pre-and post-merger will vary depending on the model, which we analyse below.
Consider two identical black holes in the merged core during the black hole merger. First, the total region of influence of the unmerged black holes will be between V = 4 3 πr 3 inf and 2V depending on separation of the black holes. Since r inf is independent of black hole mass, the total region of influence after the merger will simply be V . In all cases, E ∝Ṁ , so the amount of feedback energy available prior to the merger will be proportional toṀBH 1 +ṀBH 2 ; after the merger, we will have energy proportional toṀBH 1 +BH 2 ; if we defineṀ to be the accretion rate of one black hole prior to the merger, then for each accretion rate we studied, across the merger we have We can clearly see that the only total accretion rate in the core that remains unchanged by the merger isṀ Edd ; the total Bondi accretion rate doubles and the total viscous and drag accretion rates halves. Thus, in Models SDH, BS and WT (which use Bondi accretion), there is more feedback energy available after the merger than before, and this increased amount of energy will be distributed to a smaller region. Likewise, there is less feedback energy available in Models ONB and DQM after the merger, which will also be distributed within a smaller volume. Thus, we can argue that the galactic core becomes more energetic assuming Bondi accretion, but less energetic assuming viscous or drag accretion when two black holes merge.
In the models, Model DQM has the least stringent merger prescription, with only a separation criteria. As expected, this is the first model to merge. The most stringent Low Resolution Figure 10. Global star formation rate for the fiducial (top) and low resolution (bottom) runs. The vertical lines represent the times of first periapsis, apoapsis, second periapsis and core merger for the dynamic black hole mass models, as given in Table 4.
merging criteria is in Model BS, which is never met. The circular velocity at the more massive black hole's smoothing length is typically ten times smaller than the local sound speed, and the relative velocity between the black holes never drops below this threshold. Moreover, the chaotic motion and two-body interactions of the black holes typically keeps them further apart than the separation criteria.

Global star formation rates
During a major merger, there is typically an burst of star formation at apoapsis and again during core merger due to the infall of gas on to the galactic core (e.g. Mihos & Hernquist (1996); SDH05; DQM11). There are several characteristics that can modify the star formation rate at these epochs, including initial conditions and the inclusion of a galactic bulge and/or AGN feedback (SDH05). In Hydra, the star formation algorithm differs from those used in SDH05, BS09, ONB08 and DQM11, so we do not expect an exact reproduction of their results; however, any differences in SFRs amongst our models will be solely the result of the AGN feedback algorithm and its interaction with the star formation algorithm. The time-averaged SFR for the fiducial (top panel) and the low resolution (bottom panel) are plotted in Fig. 10. In these models, there are notable star formation bursts shortly after the beginning of the simulation, at apoapsis and again at core merger. In all cases, the initial burst is a result of the gas disc relaxing (Williamson & Thacker 2012); this initial burst also exists in simulations of an isolated galaxy that is constructed in a similar manner to the one presented here. SDH05 and DQM11 create the vertical structure of their gas disc by assuming hydrostatic equilibrium after setting their equation of state. Although cooling is not considered during this process, it means that their initial conditions and simulations are more equivalent than ours. Thus, as expected, this initial burst of star formation is missing from their results. This initial burst is resolution dependent, being higher and shorter lived in the fiducial resolution models. After this epoch in the low resolution models, the only outstanding features are post-apoapsis burst in Model SDH l , the core merger burst for Models SDH l and DQM l , and the expected drop in SFRs for all models (except Model ONB l ) shortly after core merger. The low resolution features do not necessarily have counterparts in the fiducial resolution simulations and vice versa, as a result of the stellar feedback being resolution dependent: low resolution star formation events extend over a considerably larger volume and mass than fiducial resolution events due to the larger gas particle mass and lower spatial resolution. Thus, a low resolution event has an unavoidably greater impact on its environment, and causes comparatively large regions of the galaxy to undergo a cooling delay. This explains the delayed star formation peak after apoapsis in Model SDH l , the lack of peak in Model WT l after apoapsis, and generally lower SFRs at core merger.
In Model SDH, the SFR increases by 0.9 dex at both apoapsis and core merger; in SDH05, there is no burst at apoapsis and there is an increase of 1 dex at core merger. This burst at core merger exists in all the models, except Model WT, where the high feedback rate of thermal energy prevents the gas from cooling in to the regime where star formation can commence. After the core merger burst, the SFRs drop, as in SDH05, due to the heating and total disruption of the system.
In Model WTh, the star formation burst at apoapsis peaks approximately 50 Myr earlier than in Model WT. The apoapsis burst is a result of a nuclear inflow of gas, and the lower value of hmin allows higher densities to be resolved sooner; thus inflowing gas enters the permissible density regime of star formation for Model WTh before Model WT. For the same resolution argument, Model WTh also has a small burst of star formation (an increase of 0.6 dex) at core merger. This difference in the SFR highlights one of the well known drawbacks of resolution dependency in the star formation algorithm, although developing a model that is independent of resolution awaits a full understanding of star formation in a global context.
The SFR in Model ONB remains relatively smooth until a small increase of 0.6 dex at core merger. Even after core merger, the SFR remains relatively unchanged due high core densities and moderate core temperatures.
In DQM11, there is a major star formation burst shortly after first periapsis and again starting at second periapsis (see their fig. 1, middle panel); Model DQM reproduced the burst at second periapsis. The hot ring around the void means that the star formation must occur in the galactic disc, which results in the lack of peak after first periapsis. At second periapsis, the merger is efficient enough to cre-ate a galactic core, causing the core merger star formation burst. After this final burst, the SFR drops to small values, agreeing with DQM11.

FINAL STATES
We evolved our galaxy merger simulations for 1.5 Gyr, and the characteristics of the remnants are presented here. The exception is the fiducial resolution Model ONB, which was evolved for 1.1 Gyr. From our analysis of Model ONB l , we do not expect any significant evolution between 1.1 and 1.5 Gyr. Therefore, for the final remnant of Model ONB, we present a combination of data from 1.1 Gyr, data extrapolated to 1.5 Gyr, and data from Model ONB l ; we will clearly state where the data comes from in each instance.
The final remnant in all models is a reformed gas and stellar disc surrounded by a fragmented hot gas halo. The radial profiles for gas temperature, gas density and stellar column density (averaged over all azimuthal and polar angles) for each remnant are similar, as seen in Fig. 11; the values typically span 0.5 dex at any given radius. In Figs. 12 and 13, we have plotted the final gas temperature, gas column density and stellar column density for frames of 100 and 20 kpc per side, respectively.

Gas properties
By the end of the simulation, the gas that was originally in a galactic disc is heavily depleted while the majority of the gas that was originally in the halo remains; the fraction of remaining gas for each fiducial resolution model is given in Table 7. The majority of the depleted gas has been converted in to stars as opposed to being accreted on to a black hole. The gas that is ejected from the disc has temperatures and pressures a few orders of magnitude greater than the halo gas, thus it can easily escape from the system. Therefore, the gas halo does not play a roll in recycling the ejected gas in to stars.
The amount of substructure varies considerably amongst remnants. Model ONB's remnant is a well formed disc in a smooth halo, while Models BS and WT have a loosely reformed disc in a fragmented halo. Shortly after core merger, Model BS undergoes an outburst event which expels most of the gas, and by 1.5 Gyr, this gas is starting to recondense on to the core. In Model WT, the remaining gas cools and fragments, and then begins to recondense on to the core. Thus, in both models there are similar remnants at 1.5 Gyr but different histories prior to this time. In all cases, a small disc begins to reform; the surface density profile for the reformed fiducial resolution discs is plotted in Fig. 14. The scale lengths of the discs embedded within the stellar remnant are approximately 0.5 kpc, compared to the initial value of 2.46 kpc. There are voids in the centre of the discs in Models SDH, BS and WT, thus the peak surface density occurs at 0.7-1.5 kpc from the centre. The voids persist from the residual angular moment of the gas and from the weak (but remaining) AGN feedback.

Stellar remnant
For each remnant, we calculate the stellar velocity dispersion, σ, around the black hole. To mimic observations, we calculate this using σ 2 = Re 0 σ 2 los (r)I(r)r dr where σ los is the line of sight velocity dispersion, I(r) is the projected 2D stellar density profile, and Re is the half-light (mass) radius; this is similar to what is done numerically in DQM11 and observationally in Gültekin et al. (2009). Since the final stellar configuration is highly triaxial (see the two right-most columns in Fig. 13), the σ we present in Fig. 15 is averaged over 1000 random lines of sight; we have also included the full range of σ's calculated, as well as the preferentially chosen lines of sight along the ±x−, ±y− and ±z−directions.
For the fiducial (low) resolution models, the velocity dispersion spans a range of ∼75 ± 15 km s −1 (∼90 ± 15 km s −1 ), and is latitude dependent. Since the remnant reforms a disc, there is a large σ if the line of sight lies in the plane of the disc, and σ decreases as the elevation increases. Except for Models BS and ONB, the average σ lies is within the one-sigma scatter of the observational relation.
In Model BS, the black holes never merge and do not fall within one-sigma of the observed relation; the black holes in Model BSw merge, but still do not fall within one-sigma of the observed relation. This indicates that the difference from observational expectations is not a result of the lack of  Table 7. Ratios of final to initial global mass components for each of our fiducial resolution models; black hole masses are dynamic masses (except for Model DQM) since gas is removed in discrete amounts of mgas. Gas labelled as disc gas or halo gas is based upon its initial position. We define Macc as the total mass of gas (dynamically) accreted by black holes and ∆Mgas = M gas,i − M gas,f . Depleted gas that is not accreted on to a black hole has been converted into stars. The values for Models ONB and ONBc are extrapolated to 1.5 Gyr based upon data from the final output and their low resolution counterparts. Specifically, we assume an average SFR of 1.0 M yr −1 (1.6 M yr −1 ) andṀ BH = 5.6 × 10 −3 M yr −1 (6.0 × 10 −8 M yr −1 ) for the remainder of the simulation for Model ONB (ONBc).  Figure 15. The M BH -σ relation for our simulations, along with the observed relation (red) and the one-sigma scatter (dashed) from Gültekin et al. (2009). For our five primary models at each resolution, the dot represents the average σ of 1000 random lines of sight, the horizontal bars represent the range of all calculated σ's, and the remaining three symbols on the horizontal line represent σ taken along the ±x−, ±y− and ±z−lines of sight. For our variant models, we have only plotted the average σ of 1000 random lines of sight (solid triangles). We have explicitly labelled which points/set of points corresponds to which models. The black dots in the red rectangles represent the black holes in Model BS; the lower two points represent the actual black holes, and the upper point is calculated assuming one black hole with M BH = M BH 1 + M BH 2 at the centre of mass of the black hole system. The solid triangle in the red rectangle is the result for Model BSw. The points in the green box are our initial relationships. The data for Models ONB and ONBc are taken at 1.1 and 1.25 Gyr, respectively. merger for the model parameters we have implemented. In these models, between apoapsis and second periapsis, the accretion rates are lower than the other models, which hinders black hole growth. Also, randomly depositing feedback energy around the black hole, rather than returning it isotropically or kernel-weighted, could increase the velocity dispersion of the stars by increasing the anisotropy of the gas.
In Model ONB, the black hole mass is 9.5 times lower than predicted by the MBH-σ relation, assuming the stellar velocity dispersion remains constant; a mass 3.4 times larger would bring the point within the one-sigma scatter of the relation. Since the feedback energy is returned far from the galactic core, gravity is the only influence the black hole has on the nearby stars. Thus, the self-regulating feedback mechanism is never fully realised, which could explain the lack of agreement with the observed relation. This is also true for Model ONBc, but to a more extreme extent since the low accretion rate means that the black hole only negligibly grows and cannot couple to the stellar system. ) r (kpc) Figure 11. Radial profiles of our remnants, averaged over all azimuthal and polar angles. From top to bottom, we show gas temperature, gas density, and stellar density. The profile for Model ONB is taken at 1.1 Gyr.

CONCLUSIONS
This paper details a study of five different AGN feedback algorithms, which were all run using the numerical code Hydra, started from the same initial conditions, and implemented the same star formation algorithm; two resolutions were tested for each model. The AGN feedback algorithms used in Models SDH, BS, ONB and DQM were previously found in the literature and the algorithm used in Model WT was developed specifically for this study. We also ran four additional models, each of which was a slight variation of one of the primary models to study a specific aspect of the model.
We have performed a full analysis of all models, but only presented fiducial resolution results of the five primary models and selected other results. Although the low resolution models run considerably faster, their results are noisy, and some of the physics is essentially lost. However, the low resolution studies are a good test of the varying feedback models, and provide a first-order comparison between the models. Furthermore, the continued implementation of cosmological simulations mean these resolutions will inevitably be used in the future.
By comparing these models using the same initial conditions, star formation algorithm and numerical code, any differences obtained are a direct result of the AGN feedback. Although we have tried to isolate each component of the AGN feedback algorithm, we accept that they are all intimately intertwined, and the effects cannot necessarily be disentangled from one another. We also accept that every simulation we run involves free parameters; although we have -as best as possible -matched them to their source simulations from the literature, it is possible that we can modify these parameters such that all simulations will return similar final states. Since the goal of this paper was to compare different feedback models and not to analyse any given model in great detail, we defer to the source literature for discussions on the sensitivity of the results to the free parameters. We thus refrain from making strong statements about the 'accuracy' of given models. Our principal conclusions of the five key components are as follows: (i) The black hole advection algorithm plays a key role since small displacements can cause great changes in the accretion rate. Coupling the black hole to a gas particle can allow for oscillations or chaotic motion if a void is formed around the black hole, or if the nearby gas does not meet the velocity criteria. Using a tracer mass yields smooth move-ment of the black hole, but the evolution time of the merger is decreased and disc morphology is altered in non-trivial ways such as preventing the formation of a bar. This obviously impacts other physical properties, most notably the SFR as well. Clearly evolving the mass of the black hole in a way that does not overly impact the expected evolution of the disc is an important goal. The algorithms that limit the distance a black hole can be displaced per iteration and where the direction is based upon the local stellar or total mass distribution appear to be optimal.
(ii) We test four different accretion rates:ṀB(α = 100), MB(α ≡ α(nH)),Ṁ drag andṀvisc. Each model yields a different accretion history, with substantially different qualitative profiles and quantitative differences of up to a factor of three orders of magnitude at any given time. The total black hole mass in the remnants varies by factor of 6.9, with Model BS being the least massive, and DQM being the most massive. The three models that fall within one-sigma scatter of the MBH-σ relation have final black hole masses that differ by less than 30 per cent.
(iii) The form of feedback can drastically modify the large-and small-scale systems. Delivering the feedback to the halo gas leads to little modification of the core region, resulting in high core densities and nominal temperatures. Thermal feedback delivered to the core region can drive strong outflows. Kinetic feedback delivered to the core region efficiently creates a void around the black hole, resulting in low gas temperature and density. This feedback is persistent and efficient, keeping the core properties nearly constant for all time; only a cataclysmic event, such as a core merger, can modify these characteristics.
(iv) We tested three categories of particle accretion algorithms: stochastic-unconditional, stochastic-conditional and continual-conditional; two of the three continual-conditional algorithms randomly selected which particles were accreted but never how many or how often. The continual-conditional cases gives the best agreement between the dynamical and internal mass; here, the discrepancy is never more than 2mg. The stochastic-unconditional algorithm of Model SDH, however, contains discrepancies up to 15mg. Thus, for agreement between dynamical and internal masses, continualconditional algorithms appear optimal.
(v) In Models BS and BS l , the black holes never merge; in these models, there is considerable chaotic movement of the black holes in the remnant core, preventing both merger criteria from being simultaneously met. In the remainder of the models, the black holes merge during or shortly after core merger, as one would reasonably expect. The importance of a 'reasonable' merger time is that the amount of feedback energy available increases across the merger for Bondi accretion and decreases for drag and viscous accretion.
The models presented here do not represent an exhaustive list. For example, an accretion rate suggested by Hobbs et al. (2012) is an interpolation between Bondi accretion and a free-fall accretion rate. Power et al. (2011) suggest using an accretion disc particle; in this two-tier accretion process, gas is first accreted on to an accretion disc, and then the gas is accreted from the disc on to the black hole using a viscous time-scale. The accretion disc particle method of accretion has recently been implemented in major merger simulations similar to those presented here (Wurster & Thacker 2013). Lastly, while it is possible that free parameters in different models can be adjusted such that all yield similar final states, it is very clear that this will not produce agreement throughout the evolution. While we have a vast amount of observational data of many stages of major mergers, including final remnants, it is difficult to build these various snapshots into an obvious picture of the evolutionary processes, and the fundamental theoretical development of AGN modelling still has much to contribute in this regard. As we have shown here, given the many different models of AGN feedback, building an understanding of the evolution of mergers is an important and challenging task.