# Off the beaten path: a new approach to realistically model the orbital decay of supermassive black holes in galaxy formation simulations

M. Tremmel, F. Governato, M. Volonteri, T. R. Quinn

► To cite this version:

M. Tremmel, F. Governato, M. Volonteri, T. R. Quinn. Off the beaten path: a new approach to realistically model the orbital decay of supermassive black holes in galaxy formation simulations. Monthly Notices of the Royal Astronomical Society, 2015, 451, pp.1868-1874. 10.1093/mnras/stv1060. insu-03644882

HAL Id: insu-03644882

https://insu.hal.science/insu-03644882

Submitted on 20 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Off the beaten path: a new approach to realistically model the orbital decay of supermassive black holes in galaxy formation simulations

M. Tremmel,<fnref n="1" />* F. Governato,<fnref n="1" /> M. Volonteri<fnref n="2" />,<fnref n="3" /> and T. R. Quinn<fnref n="1" />

Accepted 2015 May 8. Received 2015 April 13; in original form 2015 January 27

## ABSTRACT

We introduce a sub-grid force correction term to better model the dynamical friction experienced by a supermassive black hole (SMBH) as it orbits within its host galaxy. This new approach accurately follows an SMBH’s orbital decay and drastically improves over commonly used ‘advection’ methods. The force correction introduced here naturally scales with the force resolution of the simulation and converges as resolution is increased. In controlled experiments, we show how the orbital decay of the SMBH closely follows analytical predictions when particle masses are significantly smaller than that of the SMBH. In a cosmological simulation of the assembly of a small galaxy, we show how our method allows for realistic black hole orbits. This approach overcomes the limitations of the advection scheme, where black holes are rapidly and artificially pushed towards the halo centre and then forced to merge, regardless of their orbits. We find that SMBHs from merging dwarf galaxies can spend significant time away from the centre of the remnant galaxy. Improving the modelling of SMBH orbital decay will help in making robust predictions of the growth, detectability and merger rates of SMBHs, especially at low galaxy masses or at high redshift.

**Key words:** methods: numerical – galaxies: evolution – galaxies: kinematics and dynamics – quasars: supermassive black holes.

## 1 INTRODUCTION

Supermassive black holes (SMBHs) represent a crucial, though still poorly understood, aspect of galaxy evolution theory. SMBHs with as much as $10^9 M_{\odot}$ in mass power luminous quasars observed at $z > 6$ (Fan et al. 2001; Mortlock et al. 2011), while in the local Universe, SMBHs are ubiquitous in both massive galaxies (e.g. Gehren et al. 1984; Kormendy & Richstone 1995; Kormendy & Ho 2013) and small, bulgeless disc galaxies (Filippenko & Ho 2003; Shields et al. 2008) as well as dwarfs (Reines et al. 2011; Reines & Deller 2012; Reines, Greene & Geha 2013; Moran et al. 2014). Empirical scaling relations between the mass of SMBHs and the stellar mass and velocity dispersion of their host galaxies are indicative of co-evolution (Häring & Rix 2004; Gültekin & et al. 2009; Kormendy & Ho 2013; Schramm & Silverman 2013). Powerful outflows resulting from feedback from an active nucleus have been observed (Feruglio et al. 2010; Alatalo et al. 2011) and are thought to play a critical role in the quenching of star formation in massive galaxies and the exponential cutoff of the galaxy luminosity function at high masses (Springel, Di Matteo & Hernquist 2005b; Bower et al. 2006; Croton 2009; Teyssier et al. 2011).

Numerical simulations have proven to be a critical tool for understanding the formation and evolution of galaxies, and extensive work has already been done to implement SMBHs into these simulations (e.g. Di Matteo et al. 2003; Hopkins et al. 2005; Sijacki, Springel & Haehnelt 2009). Due to the scale of these simulations, which often include a full cosmological environment, the resolution is necessarily limited. Accretion on to SMBHs and the following feedback processes are hence implemented via sub-grid prescriptions, under the broad assumption that conditions at the smallest resolved scale drive the BH evolution at the scale of just a few parts. One major obstacle is that it becomes necessary to accurately follow the dynamics of a single object (the SMBH), something that cosmological simulations with force resolution larger than a few pc are inherently not well equipped to do. These dynamics can have important consequences for the growth of SMBHs, particularly during galaxy mergers. Detailed simulations of galaxy mergers show that the dynamics of SMBHs are intimately connected with accretion and vice versa (e.g. Callegari et al. 2011).

In this paper, we describe our approach to take into account unresolved dynamical friction (DF) in the orbital evolution of SMBHs and show that our prescription is a promising and realistic alternative to advection that can easily be implemented in existing codes. We show that this approach, when applied to simulations with the typical resolution of modern cosmological simulations (~0.1–0.5 kpc),

<footnote n="1">¹Astronomy Department, University of Washington, PO Box 351580, Seattle, WA, 98195-1580, USA</footnote>

<footnote n="2">²Institut d'Astrophysique, UMR 7095 CNRS, Université Pierre et Marie Curie, 98 bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="3">³Department of Astronomy, University of Michigan, Ann Arbor, MI, 48109-1107, USA</footnote>

<footnote>* E-mail: m.tremmel6@gmail.com</footnote>
leads to more realistic dynamics that match well with analytic approximations for DF time-scales. In Section 2, we summarize some current methods for correcting SMBH dynamics in cosmological simulations and why a new method is needed. In Section 3 we present our method. In Section 4 we present our results, comparing our method with a popular advection prescription in both an isolated dark matter (DM) halo and a fully cosmological zoomed-in simulation of a dwarf galaxy. In Section 5, we summarize and discuss our results.

## 2 THE NEED FOR A NEW MODEL

DF, the force exerted by the gravitational wake caused by a massive object moving in an extended medium (Chandrasekhar 1943; Binney & Tremaine 2008), causes the orbits of SMBHs to decay towards the centre of massive galaxies (Governato, Colpi & Maraschi 1994; Kazantzidis et al. 2005). In cosmological simulations, modelling DF is particularly challenging, as the particle representing an SMBH is often only a few times the mass of the background particles. In this case, SMBH orbits can become dynamically heated due to the limited mass resolution and the resulting spurious collisionality from a noisy gravitational potential (Hernquist & Barnes 1990). This numerical heating can cause the black holes (BHs) to gain/lose energy and be unrealistically perturbed away from the centre of its host halo. In order to lessen collisionality and improve performance, N-body simulations employ a gravitational softening length, a characteristic scale below which the gravitational force between two particles becomes damped. This mechanism, while necessary numerically, hinders DF by preventing the close interactions that are necessary to form a wake in the vicinity of the SMBH.

To account for the unrealistic dynamics of SMBHs in cosmological simulations, many groups employ artificial advection schemes that change the motion of SMBHs. Each method comes with its own drawbacks. For example, placing the SMBH at the position of the lowest potential gas particle around it (e.g. Springel, Di Matteo & Hernquist 2005a; Booth & Schaye 2009) causes chaotic motions, especially when relative velocity constraints are put on the gas particles (Wurster & Thacker 2013). Utilizing a large 'tracer mass' with which the SMBH gravitationally interacts with its surroundings (Debuhr, Quataert & Ma 2011) affects the morphology in the galactic centre as well as the accretion history of the SMBH (Wurster & Thacker 2013). Pushing the SMBH a certain distance in the direction of either the local stellar density gradient (Okamoto, Gao & Theuns 2008) or the local centre of mass (COM; Wurster & Thacker 2013) can avoid chaotic motions but adds an additional free parameter, the distance the SMBH is pushed each step. Even disregarding their individual drawbacks, none of these solutions are ideal, as they all require the explicit assumption that SMBHs rapidly decay into and remain stable at the centre of their host galaxies, which is often not accurate (Bellovary et al. 2010).

In strongly interacting systems, where SMBHs are thought to go through much of their growth and possibly their formation (Mayer et al. 2007), the inner regions of a galaxy are likely to experience strong potential fluctuations, affecting the dynamics of the SMBH. This could be especially evident in dwarf galaxies, where BHs exist (Reines et al. 2013; Moran et al. 2014) within a shallow potential with a cored density profile (Oh et al. 2011; Teyssier et al. 2013), making perturbations to their orbits more likely. Additionally, SMBHs accreted during both major and minor mergers, especially dry mergers (Kazantzidis et al. 2005), do not immediately sink to the centre of the galaxy, creating a population of small 'wandering' SMBHs (Islam, Taylor & Silk 2003; Volonteri

<table><caption>**Table 1.** Information about the resolution of each isolated DM halo test simulation.</caption><thead><tr><th>Name</th><th>Number of particles</th><th>DM mass<sup>a</sup> (M<sub>☉</sub>)</th><th>Softening length, ε<sub>g</sub> (pc)<sup>b</sup></th></tr></thead><tbody><tr><td>Low Res</td><td>1.05 × 10<sup>6</sup></td><td>9.78 × 10<sup>5</sup></td><td>311</td></tr><tr><td>Oversampled</td><td>8.39 × 10<sup>6</sup></td><td>1.22 × 10<sup>5</sup></td><td>311</td></tr><tr><td>High Res</td><td>6.71 × 10<sup>7</sup></td><td>1.53 × 10<sup>4</sup></td><td>77</td></tr><tr><td>Small Soft</td><td>6.71 × 10<sup>7</sup></td><td>1.53 × 10<sup>4</sup></td><td>10</td></tr></tbody></table>

*Notes.* <sup>a</sup>Mass of the DM particles in the simulation.

<sup>b</sup>Spline kernel gravitational softening.

& Perna 2005; Bellovary et al. 2010). Advection, by assuming that orbital decay time-scales are always short, could then artificially increase the mass of central SMBHs, through inflated merger and accretion rates.

More advanced approaches have been explored. For example, Dubois et al. (2013) include a sub-grid prescription for gas DF acting on the SMBH. While this is a promising method, it requires assumptions about the multiphase nature and equation of state of the gas far below the resolution limit of any cosmological simulation (~1–5 pc). Lupi, Haardt & Dotti (2015) increase the resolution around SMBHs to attain very accurate dynamics, but their method is applicable only for very high resolution simulations meant to closely follow BH–BH mergers, not cosmological simulations.

The method we propose in this paper is a simple solution to correcting SMBH dynamics in cosmological simulations. Our approach is to estimate the unresolved DF felt by the SMBH and apply the appropriate force to the SMBH particle. This is a significant improvement over artificial advection, as it makes no assumption about where the SMBHs should be located in a galaxy at a given time and it has no explicit effect on the SMBH surroundings. In addition, it requires minimal assumptions about the state of the simulation below the resolution limit and it naturally converges with increasing resolution.

## 3 METHODOLOGY

### 3.1 The test simulations

To test the dynamics of BHs in a realistic DM halo, we generate initial conditions for a collapsing overdensity following the procedure of Evrard (1988).<fnref n="1" /> We begin the simulation by approximating the state of a halo beginning to collapse at high redshift. We then allow the halo to collapse before implanting a BH of mass $10^6$ M$_{☉}$. We run a suite of simulations at different mass and spatial resolutions (see Table 1) with a massive ‘BH’ particle initially (1) at the centre while the halo is still actively collapsing or (2) on an eccentric orbit after the halo has mostly relaxed. The resolutions of the ‘Low Res’ and ‘Oversampled’ runs are comparable to current uniform volume simulations (Kereš et al. 2012; Dubois, Volonteri & Silk 2014a; Sijacki et al. 2014; Schaye et al. 2015), while the ‘High Res’ simulation is representative of high-resolution ‘zoomed-in’ simulations (Governato et al. 2012; Christensen et al. 2014; Hopkins et al. 2014). The force softening lengths adopted are typical of cosmological runs, while in ‘Small Soft’, a variant of the ‘High Res’

<footnote n="1">¹ Created using the package *ICING* created by M. Tremmel. The code is publicly available and can be downloaded at https://github.com/mtremmel/ICInG.git</footnote>
simulations, the force resolution is only 10 pc, typical of simulations of isolated binary mergers (Capelo et al. 2015) or very high resolution cosmological zoomed-in simulations (e.g. Dubois et al. 2014b).

We verified that the halo formed in the collapse has a density profile typical of cold dark matter haloes (Navarro, Frenk & White 1996). The halo we use here has (after a time consistent with $z \sim 0$) $M_{vir} \sim 2 \times 10^{11}$ M☉, $R_{vir} \sim 115$ kpc and a concentration $c \sim 4.5$. An analytic expression of the approximate time-scale for a rigid object to sink to the centre of a DM halo was calculated by Taffoni et al. (2003). Given these conditions and the mass of the test BH, the estimated DF time-scale ($\tau_{DF}$) for an eccentric ($v = 0.1$ $v_{circ}$) orbit at an initial distance of 2 kpc from the halo centre is approximately 1.8 Gyr. This initial orbit was chosen as typical of an SMBH after its parent satellite has been tidally disrupted. It has a non-negligible $\tau_{DF}$, but still much less than a Hubble time.

## 3.2 The DF prescription

DF occurs due to both large- (Colpi, Mayer & Governato 1999) and small-scale perturbations to the BH's surroundings. We consider perturbations on scales larger than the gravitational softening length, $\epsilon_g$, to be well resolved. On these scales the potential should be smooth, so long as enough particles are used so that dynamical heating is minimized.

The acceleration a black hole of mass $M_{BH}$ feels from DF due to particles of mass $m_a \ll M_{BH}$ is given by the following formula:

$$a_{DF} = -4\pi G^2 M_{BH} m_a \ln \Lambda \int d^3 \mathbf{v}_a f(\mathbf{v}_a) \frac{\mathbf{v}_{BH} - \mathbf{v}_a}{|\mathbf{v}_{BH} - \mathbf{v}_a|^3}, \quad (1)$$

where $f(\mathbf{v}_a)$ is the velocity distribution function, $\ln \Lambda$ is the Coulomb logarithm, $\mathbf{v}_a$ is the velocity of the surrounding background objects relative to the local COM velocity and $\mathbf{v}_{BH}$ is the relative velocity of the black hole. We assume that within $\epsilon_g$ from the BH the velocity distribution is isotropic, giving Chandrasekhar's dynamical friction formula (Chandrasekhar 1943):

$$a_{DF} = -4\pi G^2 M m_a \ln \Lambda \frac{\mathbf{v}_{BH}}{v_{BH}^3} \int_0^{v_{BH}} dv_a v_a^2 f(v_a). \quad (2)$$

This can be further simplified by substituting the integral for $\rho(< \mathbf{v}_{BH})$, which is the density of particles moving slower than the black hole:

$$a_{DF} = -4\pi G^2 M \rho(< \mathbf{v}_{BH}) \ln \Lambda \frac{\mathbf{v}_{BH}}{v_{BH}^3}. \quad (3)$$

The Coulomb logarithm depends on the maximum and minimum impact parameters, $b_{max}$ and $b_{min}$, such that $\ln \Lambda \sim \ln(b_{max}/b_{min})$. Because DF is well resolved at scales greater than the softening length, we set $b_{max} = \epsilon_g$ to avoid double counting frictional forces that are already occurring. For the minimum impact parameter, we take it to be the minimum $90^{\circ}$ deflection radius, with a lower limit set to the Schwarzschild radius, $R_{Sch}$:

$$b_{min} = \max(b_{90}, R_{Sch}); b_{90} = \frac{GM_{BH}}{v_{BH}^2}. \quad (4)$$

For the calculation, we use 64 collisionless particles (i.e. DM and star particles, if present) closest to the BH. We calculate the velocity of each particle relative to the COM velocity of those 64 particles. We verified that our results do not depend strongly on the number of neighbours used, although using too few particles could result in numerical noise in the calculation of this force. Since we

are explicitly assuming that the velocity distribution is isotropic, the following must be true:

$$\rho(< v_{BH}) = \frac{M(< v_{BH})}{M_{\text{total}}} \rho, \quad (5)$$

where $\rho$, the total density around the BH, is calculated by smoothing over the chosen 64 particles, i.e. $\rho = \sum_i^{64} m_i W(\mathbf{r}_{BH} - \mathbf{r}_i, h)$. $M(< v_{BH})$ is the total mass of the chosen particles that are moving slower than the BH relative to the local COM velocity and $M_{\text{total}}$ is the summed mass of all 64 particles.

The resulting acceleration (from equation 3) is added to the BH's current acceleration, to be integrated the following time step. As the spatial resolution or BH mass increases (or the velocity of the BH decreases), $b_{min}$ will become greater than $b_{max}$, in which case we claim that DF is being fully resolved and therefore the correction is not needed.

This method is not accounting for DF from gas, which can have important effects for supersonic BHs in regions where gas density dominates stars and DM (Ostriker 1999; Chapon, Mayer & Teissier 2013). This should not occur often on the scales relevant in these simulations. The centres of larger galaxies are dominated by stars, and smaller galaxies have significant star and DM fractions within the central regions (Oh et al. 2008, 2011), so this effect will only be a minor correction in most cases. DF may be overestimated within resonant DM cores where DF can become much less efficient (Read et al. 2006). This effect is secondary and mainly important when the gravitational softening length is appreciable compared to the size of the core structure. Often the orbital differences should be smaller than or similar to the resolution limit and therefore unimportant. Additionally, interactions with clumpy gas have been shown to significantly increase the time-scale for the orbital decay of SMBH binaries below ~100 pc (Fiacconi et al. 2013; Roškar et al. 2015). However, this effect would not be well resolved by even high-resolution 'zoomed-in' cosmological simulations.

# 4 RESULTS

## 4.1 Isolated DM halo

We find that, for SMBHs placed in the centre of a collapsing halo, only the 'Low Res' simulation experiences significant dynamical heating, causing the BH to be unrealistically perturbed away far from the halo centre. In the rest of the runs, the BH remains within one softening length from the halo centre for 6 Gyr and shows no sign of heating. This is due to the higher mass resolution present in those runs. Simulations including either our DF correction or an advection correction (see below) have the same results, showing little difference from simulations with no correction at all for this scenario.

Fig. 1 shows the orbital evolution of a BH initially on an eccentric ($v=0.1$ $v_{circ}$) orbit with apocentre of 2 kpc, placed in the halo after it has finished most of its collapse. The centre of the halo is defined at each step using the shrinking spheres method (Power et al. 2003). We verified that the density maximum and potential minimum coincide within much less than the force resolution. The vertical line represents the DF time-scale for the orbit derived from the analytic model of Taffoni et al. (2003), and the horizontal lines represent $\epsilon_g$ and $2\epsilon_g$. Without the DF correction, only the 'Small Soft' model, with 10 pc spatial resolution and DM particle mass almost 100 times smaller than the BH, is able to show substantial orbital decay within $2\tau_{DF}$. Implementing our DF correction results in a noticeable improvement for the orbital decay, even at the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Effects of the DF correction: distance of a test $10^6$ M$_\odot$ BH from the halo centre as a function of time at different resolutions. Dash-dotted and dashed lines indicate $\epsilon_g$ and $2\epsilon_g$, respectively, for both the ‘Low Res’ and ‘Oversampled’ models. The BH starts on an eccentric ($v=0.1 v_{\text{circ}}$) orbit with apocentre of 2 kpc. The vertical solid line represents the analytically derived time-scale for orbital decay. When the DF correction is applied, marked improvement is seen for all models except ‘Low Res’, which experiences too much dynamical heating due to the lower mass resolution. The orbits of ‘High Res’ and ‘Small Soft’ are very nearly the same once the correction is implemented, indicating numerical convergence when the DM particle mass is $\sim 10^4$ M$_\odot$.</figcaption></figure>

relatively modest resolution of the ‘Oversampled’ model, where it falls to within $2\epsilon_g$ of the centre before $2\tau_{\text{DF}}$. At higher resolution, the dynamics converge to closely match with the analytical approximation. Note that even for our highest resolution simulation, ‘Small Soft’, the DF correction causes the SMBH to sink almost 1 Gyr sooner. These are very encouraging results, as they indicate that this correction results in realistic BH orbital evolution even at resolutions attainable in large-volume simulations and it has important consequences even at the highest resolutions tested here.

Fig. 2 compares the performance of our DF prescription to that of a commonly used advection scheme used in Sijacki et al. (2007) and various other simulations. The test is done using the ‘Oversampled’ run, as this most closely resembles the resolution of a cosmological volume simulation. The BH is placed on an eccentric orbit, as in Fig. 1. The advection scheme adopted repositions the BH each time step to the position of the lowest potential particle within its 32 nearest neighbours while keeping the velocity unchanged. Not surprisingly, this results in the BH staying very close to the centre of the halo even when initially set on an off-centre orbit. The DF correction captures the more gradual orbital decay that the advection scheme completely misses. The run with no dynamical correction fails to have the BH sink within $2\epsilon_g$ even after 6 Gyr.

This is an important conclusion because these different orbital evolutions would result in drastically different accretion histories for the BH. Off-centre BHs should accrete less due to lower gas densities. Additionally, simulations that utilize advection would have BHs merge much sooner than what is predicted by their orbital decay time-scale. Our improved method should then have important implications for the growth and merger rate of SMBHs in cosmological simulations of galaxy formation.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** DF correction versus advection: results from the ‘Oversampled’ model when implementing our DF correction (blue) compared with a commonly used advection routine (red) and no correction to dynamics (cyan). Dash-dotted and dashed lines indicate $\epsilon_g$ and $2\epsilon_g$, respectively. The BH starts on an eccentric ($v = 0.1 v_{\text{circ}}$) orbit with apocentre of 2 kpc. Advection immediately pushes the off-centre BH to the centre, missing the orbital decay that our method captures well. Without any correction, the orbit decays far too slowly, remaining far ($>2\epsilon_g$) from the halo centre even after 6 Gyr.</figcaption></figure>

## 4.2 Cosmological dwarf galaxy simulation

As a first test of the dynamics of SMBHs in a fully cosmological setting, we run a high-resolution ‘zoomed-in’ simulation that results in two dwarf galaxies with masses $\sim 10^{10}$ M$_\odot$ at $z=0$. The
simulation has a resolution similar to our 'High Res' isolated halo model, with DM particle mass $1.6 \times 10^4 M_{\odot}$, gas particle mass $3.3 \times 10^3 M_{\odot}$ and gravitational softening of only 87 pc. We showed in the previous section that at this resolution our DF prescription gives results that match analytic models. We chose a dwarf galaxy for this test because SMBHs are more likely to become perturbed away from the galactic centre, given their shallow gravitational potential and actively evolving, cored DM profile (Governato et al. 2012; Pontzen & Governato 2012). This will guarantee a useful test environment for exploring the differences between our method and advection.

This test is also topical, as there is a growing sample of dwarf galaxies with detected SMBHs (Reines et al. 2013; Moran et al. 2014). Realistic numerical studies of BH formation and growth in these small galaxies, focusing on their occupation fraction and how they and their host galaxies evolve towards the correlation with the stellar velocity dispersion, would provide vital constraints on BH seed masses and early growth mechanisms (Volonteri, Lodato & Natarajan 2008; Volonteri & Gnedin 2009; Volonteri 2010).

We use the new N-body + SPH code CHANGA (Menon et al. 2015), which includes all of the physics modules previously implemented in GASOLINE (Wadsley, Stadel & Quinn 2004) such as hydrodynamics, gas cooling, a cosmic UV background, star formation and SNe feedback. The 'zoomed-in' approach preserves the large-scale tidal field while allowing us to model a small region at high resolution. This is similar to the simulation described in Governato et al. (2010) and in Shen et al. (2014). In this particular simulation, stars and DM densities dominate that of gas within the inner regions by more than a factor of 10. Examining other dwarf galaxy simulations (e.g. Governato et al. 2010), we find that gas densities can often reach similar values to DM and stars. In either case, the contribution from gas to DF is only a minor correction, at most a factor of a few, and negligible in our current example.

By $z < 1$ these systems are a good representation of real dwarf galaxies, with an extended stellar disc, no bulge, a high gas fraction and a cored DM profile. We first run the simulation until $z = 6$, when we insert five BHs of mass $5 \times 10^5 M_{\odot}$ in the centres of the five most massive haloes at the time, which all have a mass of $2 \times 10^8 M_{\odot}$ or higher. In these simulations, BHs do not accrete or produce feedback, as we are only interested in following their dynamics. From $z = 6$ we run two simulations to $z < 1$, one with our DF routine and the other with advection. By the end of the simulation, there are only two major star-forming galaxies with masses $\sim 10^{10} M_{\odot}$.

We then run the Amiga Halo Finder (AHF; Knollmann & Knebe 2009) for all the saved snapshots and calculate the centre of the main halo at each step using the shrinking spheres approach. In Fig. 3, we follow the trajectories of four BHs with respect to the centre of this halo (which originally just has the blue BH at its centre). Each colour represents a different BH. Black dots indicate a BH merger, which happens when two BHs come within $2\epsilon_g$ of one another at relative speeds low enough to be gravitationally bound. The dashed black lines indicate the gravitational softening length of 87 pc.

The BHs become perturbed as the red, cyan and blue host galaxies interact between $z = 3$ and $4$. In the advection case, the BHs are driven quickly towards the centre where they all merge, leaving only one BH (labelled as red). With the DF correction, only the cyan and blue BHs merge. After the red and blue galaxy hosts merge with DF, the blue and red BHs remain orbiting around the centre of the merger remnant. The blue BH comes back to the centre only after 4 Gyr and the red BH remains orbiting at around 1 kpc (11 times

<figure><img src="image_4.png" /><figcaption>**Figure 3.** DF correction in a cosmological dwarf simulation: the dynamics of four BHs in the cosmological zoomed-in dwarf galaxy simulation with DF (top) and advection (bottom). These are the BHs that end up in the most massive system by the end of the simulation. Each coloured line traces the distance of a BH from the centre of the most massive halo. Black dots mark merger events and the dashed lines mark the gravitational softening length of the simulation (87 pc). Which of the two BHs emerges from a merger event and which is 'eaten' is unimportant. DF is able to sustain a long-lived dual BH system (blue and red) while the advection scheme causes them to quickly merge. The green BH remains on a very wide orbit in the DF run, but is quickly and unrealistically pulled to the centre with advection.</figcaption></figure>

the force resolution) for another 4 Gyr before sinking and merging with the blue BH.

The more striking difference between the DF and the advection run involves the green BH. When the much smaller green host galaxy merges with the blue/red host at $z \sim 1$, it is initially far from the halo centre ($\sim 30$ kpc) and is quickly disrupted by the main galaxy. With DF, the BH stays on a wide orbit, never coming much closer than a few kpc from the halo centre. In the advection case, however, the green BH is quickly pushed to the centre where it merges with the central red BH (see Figs 3 and 4). This is an unrealistic result, as the DF time-scale of a $5 \times 10^5 M_{\odot}$ BH that far away from the centre of such a small galaxy would be longer than a Hubble time.

With this simulation we can clearly see how the choice of dynamical correction can affect the ability of SMBHs to become perturbed during mergers. In the DF simulation, BHs are able to remain off-centre for many Gyr while with advection they are quickly driven to the centre. Additionally, the DF simulation allows for sustained wide orbits resulting from minor mergers. Such dynamics can have an important impact on interpreting the connection between the initial occupation probability of SMBH seeds in dwarf galaxy progenitors and the observed occupation at low redshift. Methods such as the advection scheme presented here would predict a more direct connection, while the simulation with DF indicates that the nature of the mergers (i.e. mass ratio and orientation) can have an impact on which dwarf galaxies have observable SMBH activity and when
<figure><img src="image_5.png" /><figcaption>**Figure 4.** BHs in a cosmological dwarf: a snapshot of a zoomed-in cosmological simulation of a forming dwarf galaxy at $z=0.846$. The gas density integrated along the line of sight is shown with darker colours indicating higher densities. In the DF simulation, a previous merger has created a central binary BH system (red and blue, see Fig. 3). The separation of ~1 kpc is well resolved by the simulation, which has a force resolution of 87 pc. A more recent merger has set a third BH (in green) on a wide orbit (see Fig. 3). In the image, the green BH is at its closest approach to the galactic centre. In the advection simulation, all BHs are quickly pushed to the centre, where they merge, causing the simulation to miss these more realistic BH orbits.</figcaption></figure>

that activity occurs. The DF correction could also have important implications for BH merger rates, allowing them to become more decoupled from galaxy merger rates than advection simulations. This will affect predictions of gravitational wave detections as well as estimates for recoiling BHs, which can have an important effect on observability (e.g. Madau & Quataert 2004).

This simulation is a useful illustration of the variety of different, realistic BH orbits our method allows compared to commonly used advection schemes. In future work, we will explore the dynamics of BHs in a variety of different merger events within a cosmological context.

# 5 DISCUSSION AND SUMMARY

We have introduced a sub-grid force correction term for SMBH motion based on DF. This correction allows us to better model the orbital decay of SMBHs in numerical simulations. We have shown using controlled experiments of isolated DM haloes that this addition matches analytic predictions of the orbital decay in DM haloes with resolutions attainable by large-volume cosmological simulations. We have also demonstrated that our prescription naturally converges with resolution.

This method is a significant improvement over existing ‘advection’ methods that force a short orbital decay time-scale regardless of the dynamical state of the system. When applied to a cosmological dwarf galaxy simulation, our method results in noticeably different BH dynamics compared with the advection scheme. In particular, our prescription:

(i) models the perturbation and gradual orbital decay of a central BH during a galaxy merger,

(ii) allows for long-lived dual BH systems with close (<1 kpc) orbits,

(iii) maintains a stable central BH when appropriate and

(iv) allows for sustained wide (>5 kpc) orbit BHs.

Correctly modelling the rich orbital dynamics of a BH within its host galaxy can have important consequences for its accretion history, duty cycle and observability that were previously neglected

in simplified ‘advection’ schemes. The dynamically complex and more realistic orbits allowed by our method will have crucial implications for the early growth of SMBHs, which takes place in small, rapidly growing galaxies at high redshift (Aykutalp et al. 2014). Additionally, understanding the relative importance of different accretion mechanisms throughout an SMBH’s lifetime requires the ability to accurately model its dynamics during all phases of galaxy evolution, including merger events. Implementation of DF routines such as the one presented here will improve the ability of cosmological simulations to accurately model SMBH accretion, growth and energy deposition in the IGM and, therefore, increase the ability of simulations to interpret and predict observational results. The implementation of this approach in our future cosmological volume and zoomed-in simulations represents an exciting chance to realistically study the growth and merger rate of SMBHs across cosmic time.

# ACKNOWLEDGEMENTS

FG and TQ were funded by NSF grant AST-0908499. FG acknowledges support from NSF grant AST-0607819 and NASA ATP NNX08AG84G. MV acknowledges support from SAO award TM1-12007X, from NSF award AST 1107675 and from a Marie Curie FP7-Reintegration-Grant within the 7th European Community Framework Programme (PCIG10-GA-2011-303609). Some results were obtained using the analysis software PYNBODY (Pontzen et al. 2013). CHANGA was developed with support from National Science Foundation ITR grant PHY-0205413 to the University of Washington and NSF ITR grant NSF-0205611 to the University of Illinois. Simulations were run using NCSA Bluewaters and NASA Pleiades computers. Resources supporting this work were provided by NSF PRAC Award 1144357 and the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. We thank Andrew Pontzen, Aycin Aykutalp, John Wise, Jillian Bellovary, Alyson Brooks, and Romain Teyssier for stimulating discussions and a careful reading of the manuscript.
REFERENCES

Alatalo K. et al., 2011, ApJ, 735, 88

Aykutalp A., Wise J. H., Spaans M., Meijerink R., 2014, ApJ, 797, 139

Bellovary J. M., Governato F., Quinn T. R., Wadsley J., Shen S., Volonteri M., 2010, ApJ, 721, L148

Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn. Princeton Univ. Press, Princeton, NJ

Booth C. M., Schaye J., 2009, MNRAS, 398, 53

Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

Callegari S., Kazantzidis S., Mayer L., Colpi M., Bellovary J. M., Quinn T., Wadsley J., 2011, ApJ, 729, 85

Capelo P. R., Volonteri M., Dotti M., Bellovary J. M., Mayer L., Governato F., 2015, MNRAS, 447, 2123

Chandrasekhar S., 1943, ApJ, 97, 255

Chapon D., Mayer L., Teyssier R., 2013, MNRAS, 429, 3114

Christensen C. R., Brooks A. M., Fisher D. B., Governato F., McCleary J., Quinn T. R., Shen S., Wadsley J., 2014, MNRAS, 440, L51

Colpi M., Mayer L., Governato F., 1999, ApJ, 525, 720

Croton D. J., 2009, MNRAS, 394, 1109

Debuhr J., Quataert E., Ma C.-P., 2011, MNRAS, 412, 1341

Di Matteo T., Croft R. A. C., Springel V., Hernquist L., 2003, ApJ, 593, 56

Dubois Y., Gavazzi R., Peirani S., Silk J., 2013, MNRAS, 433, 3297

Dubois Y., Volonteri M., Silk J., 2014a, MNRAS, 440, 1590

Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., 2014b, MNRAS, 440, 2333

Evrard A. E., 1988, MNRAS, 235, 911

Fan X. et al., 2001, AJ, 122, 2833

Feruglio C., Maiolino R., Piconcelli E., Menci N., Aussel H., Lamastra A., Fiore F., 2010, A&A, 518, L155

Fiacconi D., Mayer L., Roškar R., Colpi M., 2013, ApJ, 777, L14

Filippenko A. V., Ho L. C., 2003, ApJ, 588, L13

Gehren T., Fried J., Wehinger P. A., Wyckoff S., 1984, ApJ, 278, 11

Governato F., Colpi M., Maraschi L., 1994, MNRAS, 271, 317

Governato F. et al., 2010, Nature, 463, 203

Governato F. et al., 2012, MNRAS, 422, 1231

Gültekin K. et al., 2009, ApJ, 698, 198

Häring N., Rix H.-W., 2004, ApJ, 604, L89

Hernquist L., Barnes J. E., 1990, ApJ, 349, 562

Hopkins P. F., Hernquist L., Cox T. J., Di Matteo T., Martini P., Robertson B., Springel V., 2005, ApJ, 630, 705

Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

Islam R. R., Taylor J. E., Silk J., 2003, MNRAS, 340, 647

Kazantzidis S. et al., 2005, ApJ, 623, L67

Kereš D., Vogelsberger M., Sijacki D., Springel V., Hernquist L., 2012, MNRAS, 425, 2027

Knollmann S. R., Knebe A., 2009, ApJS, 182, 608

Kormendy J., Ho L. C., 2013, ARA&A, 51, 511

Kormendy J., Richstone D., 1995, ARA&A, 33, 581

Lupi A., Haardt F., Dotti M., 2015, MNRAS, 446, 1765

Madau P., Quataert E., 2004, ApJ, 606, L17

Mayer L., Kazantzidis S., Madau P., Colpi M., Quinn T., Wadsley J., 2007, Science, 316, 1874

Menon H., Wesolowski L., Zheng G., Jetley P., Kale L., Quinn T., Governato F., 2015, Comput. Astrophys. Cosmol., 2, 1

Moran E. C., Shahinyan K., Sugarman H. R., Vélez D. O., Eracleous M., 2014, AJ, 148, 136

Mortlock D. J. et al., 2011, Nature, 474, 616

Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563

Oh S., de Blok W. J. G., Walter F., Brinks E., Kennicutt R. C., 2008, AJ, 136, 2761

Oh S.-H., de Blok W. J. G., Brinks E., Walter F., Kennicutt R. C., Jr, 2011, AJ, 141, 193

Okamoto T., Gao L., Theuns T., 2008, MNRAS, 390, 920

Ostriker E. C., 1999, ApJ, 513, 252

Pontzen A., Governato F., 2012, MNRAS, 421, 3464

Pontzen A., Roškar R., Stinson G., Woods R., 2013, Astrophysics Source Code Library, record ascl:1305.002

Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14

Read J. I., Goerdt T., Moore B., Pontzen A. P., Stadel J., Lake G., 2006, MNRAS, 373, 1451

Reines A. E., Deller A. T., 2012, ApJ, 750, L24

Reines A. E., Sivakoff G. R., Johnson K. E., Brogan C. L., 2011, Nature, 470, 66

Reines A. E., Greene J. E., Geha M., 2013, ApJ, 775, 116

Roškar R., Fiacconi D., Mayer L., Kazantzidis S., Quinn T. R., Wadsley J., 2015, MNRAS, 449, 494

Schaye J. et al., 2015, MNRAS, 446, 521

Schramm M., Silverman J. D., 2013, ApJ, 767, 13

Shen S., Madau P., Conroy C., Governato F., Mayer L., 2014, ApJ, 792, 99

Shields J. C., Walcher C. J., Böker T., Ho L. C., Rix H.-W., van der Marel R. P., 2008, in Bureau M., Athanassoula E., Barbuy B., eds, Proc. IAU Symp. 245, Formation and Evolution of Galaxy Bulges. Cambridge Univ. Press, Cambridge, p. 259

Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877

Sijacki D., Springel V., Haehnelt M. G., 2009, MNRAS, 400, 100

Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G., Nelson D., Hernquist L., 2014, preprint (arXiv:e-prints)

Springel V., Di Matteo T., Hernquist L., 2005a, MNRAS, 361, 776

Springel V., Di Matteo T., Hernquist L., 2005b, ApJ, 620, L79

Taffoni G., Mayer L., Colpi M., Governato F., 2003, MNRAS, 341, 434

Teyssier R., Moore B., Martizzi D., Dubois Y., Mayer L., 2011, MNRAS, 414, 195

Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068

Volonteri M., 2010, A&AR, 18, 279

Volonteri M., Gnedin N. Y., 2009, ApJ, 703, 2113

Volonteri M., Perna R., 2005, MNRAS, 358, 913

Volonteri M., Lodato G., Natarajan P., 2008, MNRAS, 383, 1079

Wadsley J. W., Stadel J., Quinn T., 2004, New Astron., 9, 137

Wurster J., Thacker R. J., 2013, MNRAS, 431, 2513

This paper has been typeset from a *TEX/LATex* file prepared by the author.