Black holes are believed to be crucial in regulating star formation in massive galaxies which makes it essential to faithfully represent the physics of these objects in cosmological hydrodynamics simulations. Limited spatial and mass resolution and the associated discreteness noise make following the dynamics of black holes especially challenging. In particular dynamical friction which is responsible for driving massive black holes towards the centres of galaxies cannot be accurately modelled with softened Nbody interactions. A number of subgrid models have been proposed to mimic dynamical friction or directly include its full effects in simulations. Each of these methods has its individual benefits and shortcomings while all suffer from a common issue of being unable to represent black holes with masses below a few times the simulated dark matter particle mass. In this paper we propose a correction for unresolved dynamical friction which has been calibrated on simulations run with the code KETJU in which gravitational interactions of black holes are not softened. We demonstrate that our correction is able to sink black holes with masses greater than the dark matter particle mass at the correct rate. We show that the impact of stochasticity is significant for lowmass black holes inlineformulatexmath idTM0002 notationLaTeXMrm BH le 5 Mrm DMtexmathinlineformula and propose a correction for stochastic heating. Combined this approach is applicable to next generation cosmological hydrodynamics simulations that jointly track galaxy and black hole growth with realistic black hole orbits.
['Astrophysics - Astrophysics of Galaxies']
A calibrated model for Nbody dynamical friction acting on supermassive black holes
{'No Header': ', 1-21 (2023)', 'Anna Genina 1 ★ , Volker Springel 1 , Antti Rantala 1': '1 Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748, Garching, Germany  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': 'Black holes are believed to be crucial in regulating star formation in massive galaxies, which makes it essential to faithfully represent the physics of these objects in cosmological hydrodynamics simulations. Limited spatial and mass resolution and the associated discreteness noise make following the dynamics of black holes especially challenging. In particular, dynamical friction, which is responsible for driving massive black holes towards the centres of galaxies, cannot be accurately modelled with softened 𝑁 -body interactions. A number of subgrid models have been proposed to mimic dynamical friction or directly include its full effects in simulations. Each of these methods has its individual benefits and shortcomings, while all suffer from a common issue of being unable to represent black holes with masses below a few times the simulated dark matter particle mass. In this paper, we propose a correction for unresolved dynamical friction, which has been calibrated on simulations run with the code KETJU, in which gravitational interactions of black holes are not softened. We demonstrate that our correction is able to sink black holes with masses greater than the dark matter particle mass at the correct rate. We show that the impact of stochasticity is significant for low-mass black holes ( 𝑀 BH ≤ 5 𝑀 DM ) and propose a correction for stochastic heating. Combined, this approach is applicable to next generation cosmological hydrodynamics simulations that jointly track galaxy and black hole growth with realistic black hole orbits.  \nKey words: methods: numerical - quasars: supermassive black holes - black hole physics - stars: kinematics and dynamics', '1 INTRODUCTION': "Most galaxies are observed to harbour supermassive black holes (BHs) close to their centres (see Kormendy & Ho 2013, for a review). These BHs are not only fascinating objects in their own right, they are also recognized to play a decisive role in galaxy formation and evolution (Richstone et al. 1998). During periods of rapid growth through gas accretion, supermassive BHs can liberate enormous amounts of energy, shine as quasars, and potentially influence their host galaxies through this energy output (Silk & Rees 1998; Di Matteo et al. 2005). But even in low states of accretion, the associated radio activity may arrest cooling flows in early-type galaxies and explain why their star formation is quenched (Croton et al. 2006). Modelling massive galaxy formation is thus not possible without addressing the evolution of supermassive black holes at the same time. In addition, supermassive BHs are also of interest as sources of powerful gravitational waves, emitted when they merge during the hierarchical galaxy formation process (e.g. Begelman et al. 1980; Ebisuzaki et al. 1991; Jaffe & Backer 2003).  \nModern models of galaxy formation therefore aim to jointly study the evolution of galaxies and supermassive black holes. This is achieved either with simplified semi-analytic treatments (e.g. Somerville & Primack 1999; Kauffmann & Haehnelt 2000; Cole et al. 2000) or with hydrodynamical simulations (e.g. Weinberger  \net al. 2017), where the latter approach has been pioneered nearly two decades ago (Springel et al. 2005b; Di Matteo et al. 2005) and matured to a point where current cosmological hydrodynamical simulation studies such as Illustris (Vogelsberger et al. 2014), Eagle (Schaye et al. 2015a), MassiveBlack (Khandai et al. 2015), HorizonAGN (Dubois et al. 2016), Magneticum (Dolag et al. 2016), IllustrisTNG (Pillepich et al. 2018), or Simba (Davé et al. 2019) are able to successfully reproduce many galaxy and black hole scaling relations and properties, including the quenching transition of massive galaxies that is related to their embedded black holes (Nelson et al. 2018).  \nThese galaxy formation models involve, however, numerous coarse assumptions to make the problem at least approximately tractable. These simplifications first of all relate to the accretion and feedback processes that occur during the growth of black holes, which cannot be simulated in an ab initio fashion in today's cosmological simulations due to the associated huge dynamic range of spatial and temporal scales. Instead, the black holes are represented as sink particles for which the accretion, feedback and merger physics is prescribed in terms of a subgrid model. But even the basic tracking of the spatial coordinates of these sink particles, which corresponds to the orbits of the black holes, is fraught with substantial difficulty in simulations of whole galaxies. This is because the typical mass resolution of such simulations for the dark matter component is rarely better than ∼ 10 5 M ⊙ (TNG50; Nelson et al. 2019). In fact, it is often orders of magnitude worse than that, especially when large  \ncosmological volumes are studied, meaning that the physical masses of many supermassive black holes are smaller or comparable to the mass of dark matter simulation particles ( ∼ 10 4 -10 5 M ⊙ ). This immediately implies that fiducial BH sink particles will experience substantial numerical heating effects in N-body models due to discreteness (Ludlow et al. 2019a).  \nAs a consequence of these numerical influences, the orbits of numerically represented supermassive black holes can easily become unphysical. In particular, the process of dynamical friction, which drives massive black holes to the centres of their host galaxies (Chandrasekhar 1943; Binney 1977; Colpi et al. 1999), cannot be fully resolved in modern cosmological simulations. This arises both due to the coarseness in the background particle representation and the use of gravitational softening, which limits the effects of gravitational encounters of BHs with the surrounding matter. As a result, simulated BHs may not sink to the centres of galaxies on appropriate timescales, or get spuriously ejected from galaxy centres. Since the densest regions are normally found at the very centre of a galaxy, this will invariably also modify the growth rate of BHs, which in turn impacts the strength of feedback processes as well as their coupling efficiency to the host galaxy (Callegari et al. 2011). It will also affect how quickly BHs come together after a galaxy merger, form hard BH binaries, and eventually coalesce. Clearly, a failure to properly account for the orbits of supermassive BHs can be a significant stumbling block for a realistic treatment of the co-evolution of galaxies and black holes.  \nA number of methods have therefore been proposed in the literature, aimed at correcting numerical effects that lead to supermassive black holes wandering away from the centres of their galaxies. These numerical effects are, however, distinct from physical processes that lead to similar outcomes, for instance due to black hole merger recoils from asymmetric gravitational wave emission, or slingshot kicks in 3-body black hole interactions. The prevalence of these physical processes is not yet fully understood. Their outcomes are believed to be close to stochastic (Rawlings et al. 2023; Partmann et al. 2023), while their modeling would require extremely high temporal, spatial and mass resolution that is currently not achievable in simulations of large cosmological volumes. We therefore focus in this study exclusively on the problem of black hole sinking due to dynamical friction. In the following, we first motivate our work by outlining some of the existing techniques and comment on their successes and shortcomings.", '1.1 Black hole repositioning methods': 'Repositioning methods generally do not intend to model the orbital evolution of black holes, but rather aim to keep black holes near the minimum of the gravitational potential by updating their position to the location of the most bound particle within some search kernel around the black hole at each timestep (Springel et al. 2005a; Booth & Schaye 2009). Differences in implementations mainly arise in how the minimum-potential particle is selected, in particular the type of the particle: gas (IllustrisTNG; Weinberger et al. 2017), or all particle types (EAGLE; Schaye et al. 2015b). Additional constraints are sometimes invoked based on the velocity of the black hole particle compared to the local gas bulk or sound speed, to avoid black holes being dragged by gas particles with low potential, but high kinetic energy (Bahé et al. 2022).  \nRepositioning methods tend to sink black holes to the centre of the galactic potential too fast; however, despite this aggressive treatment, this method has been shown to suffer from other unphysical behaviours such as removal of black holes from backsplash galaxies  \n(Borrow et al. 2023), or strongly enhanced merging of black holes and their seeds at high redshift (Ma et al. 2021; Chen et al. 2022). Some more recent works (Pellissier et al. 2023) have introduced new repositioning models motivated by gradient descent methods that, to an extent, mimic the orbital decay of a black hole. The sinking time of the black hole can also be increased or decreased using a weight function that determines how aggressively the black hole is drifted in the direction of the potential gradient, and by limiting its velocity using a Courant-like criterion. While delaying the sinking time of the black holes compared to standard repositioning methods, the gradient descent models still cannot accurately reproduce the sinking rate of individual black holes.', '1.2 Compensating for unresolved dynamical friction': "A number of works based their subgrid dynamical friction prescription on the Chandrasekhar formula (Chandrasekhar 1943), describing the deceleration that a massive body experiences when traversing a 'sea' of isotropically distributed low-mass particles:  \na df = -16 𝜋 2 ln Λ 𝐺 2 𝑀 BH ∫ 𝑣 0 𝑓 ( 𝑢 ) 𝑢 2 d 𝑢 𝑣 3 v , (1)  \nwhere 𝐺 is the gravitational constant, 𝑓 ( 𝑢 ) is the distribution function of background particles (where often the Maxwellian distribution is assumed), 𝑀 BH is the black hole mass, 𝑣 is the black hole velocity and ln Λ = ln 𝑏 max 𝑏 min is the Coulomb logarithm. The minimum impact parameter, 𝑏 min , is the impact parameter for a 90-degree deflection. The maximum impact parameter, 𝑏 max , is not well defined but, in principle, describes the largest impact parameter where the gravitational scattering of a background particle still contributes to dynamical friction.  \nThe work of Tremmel et al. (2015) introduced an estimator aimed at compensating for unresolved dynamical friction below the gravitational softening scale,  \na df , T = -4 𝜋𝐺 2 𝑀 BH 𝜌 ( 𝑣 < 𝑣 BH ) ln Λ T v BH 𝑣 3 BH , (2)  \nwhere 𝜌 ( 𝑣 < 𝑣 BH ) is the density of particles moving slower than the black hole, measured in a kernel covering 64 nearest collisionless particles. The Coulomb logarithm is here defined as:  \nΛ T = 𝑏 max , unres ( 𝐺𝑀 BH )/ 𝑣 2 BH , (3)  \nwhere 𝑣 BH is the velocity of the black hole relative to the surrounding medium, and the maximum unresolved impact parameter, 𝑏 max , unres , is taken to be the gravitational softening, 𝜖 𝑔 . Such a correction has been employed in the Romulus cosmological simulations suite to study supermassive BH dynamics and populations of wandering black holes (Tremmel et al. 2017; Ricarte et al. 2019, 2021).  \nAs pointed out by Chen et al. (2022), the kernel over the nearest 64-100 particles can be noisy and the choice of 𝜖 𝑔 as the maximum unresolved impact parameter is not necessarily justified, given that the forces are softened compared to a Newtonian case beyond the scale of gravitational softening in many simulations (but dependent on how softening is implemented). It is also worth noting that, for particularly small kernel sizes, there may be no particles moving slower than the black hole. Chen et al. (2022) instead propose the formula:  \na df , C = -4 𝜋𝐺 2 𝑀 BH 𝜌 sph F ( 𝑥 ) ln Λ C v BH 𝑣 3 BH , (4)  \nwhere 𝜌 sph is the density of collisionless particles measured within  \nthe Smooth Particle Hydrodynamics (SPH) kernel employed in their simulations, and F( 𝑥 ) describes a Maxwellian local velocity distribution:  \nF( 𝑥 ) = erf ( 𝑥 ) -2 𝑥 √ 𝜋 𝑒 -𝑥 2 , 𝑥 = 𝑣 BH √ 2 𝜎 𝑣 , (5)  \nwhere 𝜎 𝑣 is the one-dimensional velocity dispersion of particles in the SPH kernel, and 𝑣 BH denotes the velocity of the black hole relative to the surrounding medium. The Coulomb logarithm is again defined as:  \nΛ C = 𝑏 max , unres ( 𝐺𝑀 BH )/ 𝑣 2 BH = 𝜂𝜖 𝑔 ( 𝐺𝑀 BH )/ 𝑣 2 BH , (6)  \nwhere 𝜂 is a multiple of the gravitational softening 𝜖 𝑔 below which dynamical friction is not resolved. In the simulations of Chen et al. (2022), this was found to be close to 𝜂 = 6.  \nThis method of Chen et al. (2022) has been shown to be less sensitive to noise (due to the use of the SPH kernel in computing average background particle properties) and allows one to take the scale below which the dynamical friction is unresolved as a free parameter of the simulation. The method has been employed in the cosmological simulation suite ASTRID (Ni et al. 2022; Weller et al. 2023).  \nThe method of Pfister et al. (2019) is similar to those described above, but takes into account the contribution of fast-moving particles to dynamical friction, which may be important in enhancing dynamical friction in cases where the density distribution of the background particles becomes shallow (Antonini & Merritt 2012; Dosopoulou & Antonini 2017), such as in central density cores (de Blok 2010). However, we also note that the 'standard' Chandrasekhar treatment is known to fail in reproducing the effect of stalling in cored density distributions (Read et al. 2006; Kaur & Sridhar 2018; Banik & van den Bosch 2021).", '1.3 Dynamical friction from 𝑁 -body particle distributions': "The work of Ma et al. (2023) points out a number of shortcomings in the application of the Chandrasekhar formula to model dynamical friction in simulations, namely the assumptions of an isotropic and homogeneous medium, ambiguities in defining the Coulomb logarithm, as well as difficulties in determining the local velocity distribution. The authors instead present a derivation of a discrete dynamical friction formula, allowing the full dynamical friction force to be computed self-consistently in 𝑁 -body simulations. Ma et al. (2023) also provide an expression for the perpendicular component of dynamical friction which arises due to the presence of a local density gradient, in addition to that parallel to the black hole velocity vector. However, the contribution of the transverse component is expected to be relatively minor (see also Just & Peñarrubia 2005).  \nUsing their method, the dynamical friction force can be computed straightforwardly with direct summation, or with tree methods by treating a distant node as a single particle. Manifest momentum conservation can also be introduced by applying an opposite force to each particle that contributes dynamical friction to the black hole - a feature that cannot be self-consistently implemented alongside models of Tremmel et al. (2015) and Chen et al. 2022 (as well as the model presented in this work, see Appendix C). Nevertheless, for relatively well-resolved cases, such an approach, in effect, double-counts the dynamical friction because a fraction of the force is already resolved in the simulation and taken into account in the force computation. Maet al. (2023) acknowledge this shortcoming and point out that this  \ncan be corrected in the future by some mass-resolution-dependent weight function.  \nFor low-resolution cases, the noisy nature of the forces also translates into a noisy dynamical friction estimate. Moreover, as the approach of Ma et al. (2023) takes into account gravitational softening in their computation of dynamical friction, it is implied that the impact parameters below the softening length scale are still unaccounted for. Ultimately, when defining a dynamical friction model, one needs to decide whether they aim to model the orbits of black holes (independent of mass) in the softened 𝑁 -body potential or in the fiducial infinitely high-resolution 𝑁 -body potential that the simulation is aiming to sample. In this work, we take the latter approach and aim for convergence in black hole sinking timescales across resolution levels.  \nThe recent work of Damiano et al. (2024) introduces a discrete estimator for unresolved dynamical friction, where each particle lying within the softening length of the black hole contributes to an approximation of the local distribution function and exerts a force on the black hole in the direction defined by the vector v BH -v m , with 𝑣 m the velocity of the particle. The maximum impact parameter in the Coulomb logarithm, as in Tremmel et al. (2015), is assumed to be the softening length, while the minimum impact parameter for each encounter is defined by the relative velocity of the two interacting particles, 𝑏 min = 𝐺 ( 𝑀 + 𝑚 ) ( 𝑣 BH -𝑣 m ) 2 . This method is thus able to treat interactions with multiple populations of background particles, which can differ in their mass and velocity distribution, consistently, without grouping them into a background 'sea' of particles as in the standard approach of Chandrasekhar (1943). However, a clear disadvantage of such a treatment is that its effectiveness relies on the presence of particles within the softening length. Issues may occur in the outskirts of halos, where not many particles are present, or even closer to their centres, if the gravitational softening is set too low for the halo mass resolution.", '1.4 Partial or no added dynamical friction': 'Finally, we note that in simulations like H/o.pc/r.pc/i.pc/z.pc/o.pc/n.pc-AGN no additional dynamical friction force due to collisionless particles is applied, instead relying on the dynamical friction already present in the simulations when the black hole mass is higher than that of the dark matter particles (Dubois et al. 2012, 2014) and an additional force due to gas drag (Ostriker 1999; Chapon et al. 2013). In /m.pc/a.pc/g.pc/n.pc/e.pc/t.pc/i.pc/c.pc/u.pc/m.pc simulations, the black holes are seeded at the location of the star particles with highest binding energies, ensuring stability of the black hole position against mergers and local gas flows, while the dynamical friction force following the Chandrasekhar formula is applied over all impact parameters out to the halo size for as long as the black holes are less massive than the dark matter particle mass (Hirschmann et al. 2014; Steinborn et al. 2015).', '1.5 This work': 'Perhaps the biggest shortcoming of the available subgrid dynamical friction models is a lack of model calibration on high-resolution simulations with unsoftened black hole interactions, in which dynamical friction is fully resolved. In this work, we achieve this through comparison to simulations run with the code KETJU. We chose to base our dynamical friction model on the previous works of Tremmel et al. (2015) and Chen et al. (2022), as these methods, in principle, only make assumptions about the effect of dynamical friction in the direct vicinity of the black hole, modified by gravitational softening,  \nTable 1. The black hole to dark matter mass ratio, the number of particles used, the number of particles within the original orbital radius of the black hole and the employed gravitational softening for our primary set of simulations carried out with GADGET-4. In bold, we show the resolution of the simulation carried out with KETJU.  \nwhile dynamical friction due to distant particles is already accounted for in the simulation (Spinnato et al. 2003). We do not consider the effect of gas drag in this work, which can be important for black hole dynamics, particularly at high redshift (Ostriker 1999; Chapon et al. 2013; Escala et al. 2005; Chen et al. 2022). We note, however, that our findings regarding the effect of gravitational softening on dynamical friction from collisionless particles in simulations may be applicable to the treatment of gas drag in the future.  \nThis paper is structured as follows. In Section 2 we introduce our numerical set-up and the simulations with KETJU. In Section 3 we evaluate how the methods of Tremmel et al. (2015) and Chen et al. (2022) behave in this set-up. In Section 4, we outline our calibration procedure to determine the optimal set of parameters of our dynamical friction model. In Section 5, we summarise the required steps for implementation of our model and apply it to a number of test cases. Finally, in Section 6, we present our conclusions.', '2 METHODS': 'In this section, we consider how well dynamical friction is resolved in 𝑁 -body simulations. To be able to quantify this, we must first define what the true rate of black hole orbital decay is in a particular set-up. For the latter, we will consider both the analytical prediction of the Chandrasekhar formula as well as the result of high resolution simulations with KETJU, the details of which we outline below.', '2.1 Simulations with GADGET-4': "Weset up a dark matter halo of mass 10 13 M ⊙ following a Hernquist density profile (Hernquist 1990), using the MAKENEWDISK code (Springel et al. 2005a). A black hole with mass 𝑀 BH = 10 9 M ⊙ is set on a circular orbit, starting at 20 kpc, corresponding to 1 / 4 of the scale radius of the Hernquist halo. We note that this set-up by no means represents a typical cosmological scenario, but we aimed to chose a test case such that i) a black hole starting off on a circular orbit sinks within the age of the Universe and ii) the sinking occurs over a sufficiently extended period of time such that we have a good probe of the black hole dynamics and the differences from predictions of high-resolution simulations with KETJU become clearly identifiable. The black hole must be massive enough to sink to the centre, but not too massive such that it causes significant motion in the centre of mass of the halo.  \nWhile keeping the halo and black hole mass fixed, we explore  \nthe black hole dynamics at varying dark matter mass resolution. We select an 'optimal' gravitational softening in each case, following the prescription of Power et al. (2003):  \n𝜖 opt = 4 𝑅 200 √ 𝑁 200 , (7)  \nwhere 𝑅 200 is the radius enclosing 200 times the critical density of the Universe (here we take 𝐻 0 = 100 km s -1 Mpc -1 ) and 𝑁 200 is the number of particles within this radius. The number of particles used, the black hole to dark matter mass ratio, and the optimal gravitational softening 1 are listed in Table 1.  \nOur primary N-body simulations are run with the code GADGET-4 (Springel et al. 2021), a parallel multi-purpose cosmological N-body code that can optionally also include gas with smoothed particle hydrodynamics. We employ its fast multipole method (FMM) of order 5 for computing gravity in this work, combined with a very conservative force accuracy setting so that the influence of residual force errors is negligible in the calculations carried out here.", '2.2 Analytical predictions': 'Our analytical prediction is based on applying dynamical friction with the Chandrasekhar formula (Eq. 1), where the distribution function is replaced with the Hernquist distribution function, 𝑓 ( 𝑟, 𝑢 ) (Chatterjee et al. 2002). Following Just et al. (2011), the minimum and maximum impact parameter in the Coulomb logarithm is defined as 𝑏 min = 𝐺𝑀 BH /( 𝑣 2 BH + 2 / 3 𝜎 2 ), where the denominator describes the typical velocity of encounters in the vicinity of the black hole. The maximum impact parameter is defined as 𝑏 max = 𝜌 /|∇ 𝜌 | , where 𝜌 is the local density and ∇ 𝜌 the density gradient. For a Hernquist cusp, 𝑏 max is equal to | 𝑟 -𝑟 0 | , the position of the black hole from the halo centre. We chose to use 𝑏 max = 𝑅 BH in the Coulomb logarithm. We found that this choice captures accurately the behaviour of the black hole in KETJU, while setting 𝑏 max = 𝜌 /|∇ 𝜌 | tends to sink the black hole substantially more slowly. This distance- and velocity dispersion-dependent choice of the Coulomb logarithm has been shown to yield a good approximation for point-like masses sinking in cuspy haloes with 𝜌 ∝ 𝑟 -1 (Hashimoto et al. 2003; Just & Peñarrubia 2005; Just et al. 2011; Arca-Sedda & Capuzzo-Dolcetta 2014). In practice, we set up a black hole particle on a circular orbit within an analytical Hernquist potential within GADGET-4, and apply an additional dynamical friction force using the above parameterisation of the Chandrasekhar formula.  \nOne issue with the Chandrasekhar treatment is that as the black hole approaches the centre of the halo, 𝑏 min can become greater than 𝑏 max (since 𝑏 max is defined as the black hole distance) and, as the velocity dispersion tends to zero at the centre of a Hernquist halos, 𝑏 min becomes very large. The effect, if left untreated, can lead to black holes getting kicked out of the halo centre when the Coulomb logarithm flips sign. We address this problem by switching offdynamicalfriction when this occurs. As a result, black holes do not fully sink to the centre, but settle on a circular orbit at a radius where 𝑏 max > 𝑏 min . We note that this is an issue that is also encountered by approaches like that of Tremmel et al. (2015), though in this case this would correspond to dynamical friction being fully resolved (the minimum impact parameter is greater than the scales over which dynamical friction is unresolved). Additionally, as simulated halos', 'No DF prescription': "<!-- image -->  \nFigure 1. Left: Distance of the black hole from the centre of the halo, as determined by the shrinking spheres algorithm (Power et al. 2003), as a function of time, for simulations with GADGET-4 that do not include a dynamical friction prescription. Blue to red lines show simulations varying from lowest to highest resolution, with the softening chosen optimally for each halo mass resolution. The black solid line is the analytical prediction using the Chandrasekhar formula (C43), and the black dashed line is the result of the run with KETJU. Right: as the left figure, but using a fixed softening of 0.44 kpc for all mass resolutions.  \n<!-- image -->  \nhave a central density core with size ∼ 𝜖 𝑔 , the black holes tend to settle on a circular orbit with an orbital radius similar to 𝜖 𝑔 .  \nWe note that although this approach takes into account the variation in the background density of the isotropic Hernquist halo, it still applies the dynamical friction only parallel to the black hole's motion. The latter is a consequence of assuming an isotropic and homogeneous particle background in the Chandrasekhar formula, such that perpendicular contributions cancel due to symmetry. This is not the case in reality for a Hernquist halo, where a density gradient exists. Nevertheless, as we show below, this assumption seems to return a correct black hole sinking timescales when compared to high-resolution simulations. Past works have also found that the contribution to dynamical friction due to density inhomogeneity is rather small (Just & Peñarrubia 2005).", '2.3 Predictions from simulations with KETJU': 'It is important to note that the Chandrasekhar formula makes a number of assumptions, in particular an isotropic and homogeneous nature of the particle distribution, which leads to dynamical friction acting opposite to the direction of motion, rather than also having a perpendicular term. The background particles are also assumed to move on hyperbolic orbits, such that, for example, no particles can become bound to the black hole, which could however effectively increase its dynamical mass. The self-gravity of the wake is also ignored, as well as any fundamental changes to the density field and the distribution function that the massive perturber would induce (Weinberg 1989; Banik & van den Bosch 2021). Despite these issues, the Chandrasekhar formula has been shown to provide a fairly accurate description of orbital decay in a number of set-ups (Velazquez & White 1999; Hashimoto et al. 2003; Just et al. 2011) and less successfully in others (van den Bosch et al. 1999; Boylan-Kolchin et al. 2008).  \nIn order to make sure that the Chandrasekhar description works for our set-up, we also run a high-resolution 𝑁 -body simulation using the KETJU code (Rantala et al. 2017; Mannerkoski et al. 2023). The initial conditions we use for this simulation are highlighted in bold in Table 1. The essence of the GADGET-4-based KETJU code is to add a regularized region around each massive black hole in the simulations in which the gravitational dynamics is accurately integrated without gravitational softening. Throughout this work we use the regularized region size of 2 . 8 𝜖 𝑔 in radius, resulting in region sizes up to a few hundreds of pc, ensuring that all interactions with the black hole in the simulations are always non-softened. The interactions between dark matter particles in the region remain softened. This allows us to capture the dynamical friction below the softening length, while also accounting for deviations from the Chandrasekhar approximation (e.g. like the assumption of hyperbolic orbits). Moreover, since a number of previous methods used to include unresolved dynamical friction in simulations rely on the Chandrasekhar description, we aim to learn about any potential shortcomings of these methods through this test.  \nIn KETJU, the integration around black holes (the MSTAR integrator, Rantala et al. 2020; Mannerkoski et al. 2023) relies on algorithmic regularization techniques (Mikkola & Tanikawa 1999; Mikkola & Merritt 2006), a minimum spanning tree coordinate system and the Gragg-Bulirsch-Stoer (GBS) extrapolation method (Gragg 1965; Bulirsch & Stoer 1966) to achieve an extremely high integration accuracy. The efficient two-fold parallelization of MSTAR allows for integrating up to ∼ 5000 particles in a single regularized region. In order to properly resolve dynamical friction, KETJU requires a mass ratio of 𝑀 BH / 𝑀 DM ≳ 100 between the black hole and the dark matter particles (Rantala et al. 2017). In this work, we use the mass ratio 𝑀 BH / 𝑀 DM = 5000. KETJU has been successfully applied for simulating massive black hole and black hole binary dynamics in both isolated galaxy mergers (e.g. Rantala et al. 2018) and cos-  \nn simulations (Mannerkoski et al. 2022). Note that in this work we employ only Newtonian gravity and do not include post-Newtonian terms in the equations of motion close to the black hole, something that is in principle available in KETJU, as the dark matter mass elements do not represent physical particles. For the user-given accuracy parameters of KETJU, we adopt a GBS tolerance of 𝜖 GBS = 10 -7 and an integration end-time tolerance of 𝜖 t = 10 -2 .', '2.4 Results for simulations without added dynamical friction': "The left panel of Fig. 1 shows the results of our numerical experiment to understand the extent to which dynamical friction is resolved in simulations of different mass and spatial resolution. The solid black line shows the analytical prediction of the Chandrasekhar formula (C43) and the dashed black line is the prediction of KETJU. As expected, an increase in resolution leads to a reduction in noise, as manifested by the oscillations in black hole peri- and apocentres. For 𝑀 BH / 𝑀 DM > 10 and 𝜖 < 2 kpc the black hole apocentre decreases with time. This highlights that even if the mass of the black hole is greater than that of the dark matter particle, sinking is not guaranteed.  \nIt can also be seen that an increase in mass and spatial resolution leads to a faster sinking of the black holes. As the mass resolution and the ratio 𝑀 BH / 𝑀 DM increase, the magnitude of the dynamical friction force becomes significantly larger than the impact of noise. A smaller gravitational softening also means a larger range of resolved impact parameters contributing to dynamical friction, thus also making the black holes sink faster. Nevertheless, even our highest resolution run with 10 8 particles does not converge to the analytical prediction, or the prediction of KETJU. As shown in Fig. 3 of Rantala et al. (2017), the gravitational softening of a GADGETtype code must be close to ∼ 1 pc to reach convergence with the KETJU results. It is also evident that the biggest change in the sinking timescale occurs when the softening is reduced to below 1 kpc, which substantially increases the magnitude of the resolved Coulomb logarithm. Another decrease in the sinking timescale of nearly 2 Gyr occurs as the softening is reduced to ∼ 100 pc, which again suggests that scales < 10 pc must be resolved in order to obtain converged predictions. Overall, it is unlikely that the correct sinking timescale can be captured for any softening that is greater than the typical 𝑏 min of the black hole in a given halo. In our case, this implies that softening ⪅ 20 pc is required, which is substantially less than the 'optimal' softening needed to suppress 2-body effects in dark matter, even in the case of 10 8 particles in the halo.  \nEvidently, gravitational softening plays a crucial role in resolving the effects of dynamical friction, and for sufficiently high mass resolution ( > 10 7 particles in the halo) its role is more important than the mass resolution itself. To confirm this, we additionally show the sinking of the black hole at each mass resolution with a fixed softening of 0.44 kpc on the right panel of Fig. 1. It can be seen that the run with 𝑀 BH / 𝑀 DM = 5 is now showing some black hole sinking (although the black hole is clearly unable to sink within a Hubble time), while for mass ratios > 5, black holes are able to sink to the centre. On the other hand, simulations with high mass resolution, which on the left panel of Fig. 1 had 𝜖 𝑔 < 0 . 44 kpc, now have delayed sinking times. We conclude that in the regime where the black hole particle mass is substantially higher than that of the dark matter the gravitational softening is more important than the mass resolution in resolving the dynamical friction force. The right panel of Fig. 1 also shows that for the case of 𝑀 BH / 𝑀 DM = 1, the black hole is kicked to higher apocentres when the gravitational softening is too small. This is because such a small softening is unable to effectively suppress the noise in the gravitational potential at this mass resolution, leading to the black  \nhole getting kicked outwards (note that this also applies to the dark matter particles).  \nA particularly interesting case is the simulation with 10 5 particles ( 𝑀 BH / 𝑀 DM = 10), which, with a smaller softening, seems to have its black hole sinking on roughly the correct timescale, unlike its counterparts with higher mass resolution. The enhanced 2-body scattering in background particles due to the small softening, combined with the particulars of the initial conditions, have changed the black hole's orbit in such a way that it happens to come close to the centre of the halo early on and thus also experiences more dynamical friction early on. Despite the apparent success in representing black hole sinking in this case, we interpret this as more of a lucky coincidence and a warning against an inadequate use of gravitational softening at a given mass resolution in simulations (i.e. noise can by chance cause black holes to sink faster, but such a behaviour is not desirable overall in our effort to obtain converged sinking timescales for black holes at different mass resolutions).", '3 PERFORMANCE OF EXISTING METHODS': 'In the following, we will evaluate how the existing methods to correct for unresolved dynamical friction perform on our test halo, namely the method of Tremmel et al. (2015) and Chen et al. (2022). We will not test the method of Pfister et al. (2019), which was developed for an adaptive mesh refinement code RAMSES (Teyssier 2002). For both the Tremmel et al. (2015) and the Chen et al. (2022) methods, we compute the relevant quantities (e.g. local density of dark matter particles, velocity of the black hole relative to the surrounding medium, velocity dispersion) at each gravity step and add the dynamical friction force from Eqs. (2) and (4) to the force computed by the gravity solver in GADGET-4. We also note that we do not boost the dynamical mass of the black hole for our lowest-resolution case, as would be done by a factor of 3 in the approach of Tremmel et al. (2015) and by a factor of 2 in Chen et al. (2022). In this section, we will only consider simulations with up to 𝑁 DM ≤ 10 7 particles, representative of relatively well-resolved haloes in large-volume cosmological simulations.', '3.1 The method of Tremmel et al. 2015': 'To estimate the local density, 𝜌 tot , Tremmel et al. (2015) find the 64 nearest neighbours of the black hole and use them to estimate the local density with a smoothing kernel 𝜌 tot = " 𝑁 𝑖 𝑚 𝑖 𝑊 (| 𝑟 -𝑟 BH | , ℎ ) . The mass of particles within the kernel that are moving slower than the black hole is counted and their density, 𝜌 ( 𝑣 < 𝑣 BH ) , is computed as 𝜌 tot 𝑀 ( 𝑣 < 𝑣 BH ) / 𝑀 tot . The velocities are taken as those relative to the centre of mass velocity of the neighbouring particles. In our application of the method, we instead compute the relevant quantities over 100 nearest particles to reduce the particle noise, and we take 𝑊 as the standard GADGET-4 smoothing kernel, with ℎ being the distance to the furthest neighbouring particle (see Eq. 108 of Springel 2010). We compute the centre-of-mass velocity and velocity dispersion as smoothed quantities, however, we have also verified that this choice does not noticeably affect the results of the simulations. As in Tremmel et al. (2015), 𝑏 max in the Coulomb logarithm is taken to be the gravitational softening of each simulation.  \nThe results of this experiment can be seen in the top left panel of Fig. 2, where the black holes, on average, sink slower than predicted by KETJU by almost 1 Gyr. One exception to this is the black hole with 𝑀 BH / 𝑀 DM = 10, which, due to noise, approaches within  \nFigure 2. Left: black hole trajectories at different halo mass resolutions after applying the Tremmel et al. (2015) dynamical friction correction, where the relevant quantities are computed using a weighted smoothing kernel containing 100 nearest dark matter particles. The maximum unresolved impact parameter is set to 𝑏 max , unres = 𝜖 𝑔 . Middle: the result of applying the Chen et al. (2022) dynamical friction correction, which assumes the Maxwellian approximation and 𝑏 max , unres = 6 𝜖 𝑔 . Right: applying the Chen et al. (2022) correction, while setting 𝑏 max , unres = 𝜖 𝑔 .  \n<!-- image -->  \nFigure 3. The estimated density of particles moving slower than the black hole, 𝜌 ( 𝑣 < 𝑣 BH ) , compared to the density derived from the Hernquist distribution function, 𝜌 true . The black line shows the prediction, the thin red line corresponds to the use of the Maxwellian approximation, as in Chen et al. (2022), and the thin blue line is the kernel-smoothed estimate of Tremmel et al. (2015). The thick lines are moving averages over 0.2 Gyr. Both densities are estimated using a kernel size of 100 particles in a simulation with 10 6 dark matter particles.  \n<!-- image -->  \n10 kpc of the centre of the halo already at ∼ 3 Gyr, leading to enhanced dynamical friction. Interestingly, the 𝑀 BH / 𝑀 DM = 1 black hole, which would show no noticeable dynamical friction without a correction, is also able to sink. We interpret this as a result of 𝑏 max being extremely large in this simulation (14 kpc), such that the total dynamical friction force is able to overcome the stochastic heating. Otherwise, the sinking time of the black holes is approximately converged between different resolution levels.', '3.2 The method of Chen et al. 2022': 'The approach of Chen et al. (2022) builds on that of Tremmel et al. (2015), where now the kernel size over which the distribution function is estimated is based on the SPH kernel, defined by neighbouring  \ngas particles 2 in cosmological simulations, rather than a fixed number of collisionless particles, which helps to reduce the noise in the density estimate. The maximum unresolved impact parameter, 𝑏 max , unres , is now defined as some factor 𝜂 times the gravitational softening, where the optimal value according to their experiments is 𝜂 = 6. It is certainly true that the exact value of the gravitational softening is not necessarily equal to the full extent of the unresolved dynamical friction scale. In many implementations of the gravitational softening, the forces are softened up to a distance of ∼ 2 . 8 × 𝜖 𝑔 , beyond which they are fully Newtonian (Springel 2010). As such, we expect that increasing the value of 𝑏 max above 𝜖 𝑔 will likely lead to a compensation of more of the unresolved dynamical friction and, in principle, black hole sinking rates which are more accurate.  \nThe result of applying the Chen et al. (2022) formula can be seen in the middle panel of Fig. 2, where the black holes now seem to sink, on average, faster than the prediction of KETJU by up to 1 Gyr, and there appears to be a slight resolution-dependent gradient in the sinking time. This may suggest that 𝑏 max = 6 𝜖 𝑔 in the Coulomb logarithm is too large of a dynamical friction correction, particularly in cases where dynamical friction is to a large part resolved.  \nTo understand the effect of the enhanced 𝑏 max in the Coulomb logarithm, we now repeat our test with the Chen et al. (2022) formula, but setting 𝑏 max , unres = 𝜖 𝑔 . This is shown in the right panel of Fig. 2. It can be seen that at the highest resolution of 𝑀 BH / 𝑀 DM = 1000, the sinking of the black hole is very similar to Tremmel et al. (2015), while lower mass-ratio black holes sink too slowly. Thus, the maximum unresolved impact parameter 𝑏 max , unres must be artificially increased by a factor of 6 to strengthen the effect of dynamical friction in these cases. Together with the shallower descent of the black holes that we see on the right panel of Fig. 2, this suggests that the Maxwellian approximation delays the sinking of the black holes, particularly when the resolution is low. Previous works have pointed out  \nthis inapplicability of the Maxwellian approximation to the centres of isotropic Hernquist haloes (Karl et al. 2015).  \nIn Fig. 3 we show the accuracy of the 𝜌 ( 𝑣 < 𝑣 BH ) estimates using the Maxwellian approximation (red) and using the kernel-weighted density of particles moving slower than the black hole (blue). Both are applied at each time step to the black hole with 𝑀 BH / 𝑀 DM = 100. We can see that the Maxwellian approximation is substantially more smooth than the approach of Tremmel et al. (2015) and, on average, slightly underestimates the density of slow-moving particles. This underestimation is significant when the black hole descends to below 10 kpc (> 3 Gyr).  \nWe additionally note that setting 𝑏 max , unres = 6 𝜖 𝑔 appears to be more effective in keeping the black holes in the centre of the halo (middle panel of Fig. 2, compared to the left and right panels). This is because the condition 𝑏 min > 𝑏 max , unres sets the dynamical friction correction to zero. Thus, the higher 𝑏 max , unres is, the more do the black holes move towards the local dark matter density centre. Another advantage of using the Maxwellian approximation is that it ensures 𝜌 ( 𝑣 < 𝑣 BH ) > 0 at all times, such that there is always a dynamical friction force, provided 𝑏 min < 𝑏 max , unres . We note however, that this may not be beneficial in shallow density profiles or dark matter cores, where no background particles are moving slower than the black hole, resulting in delayed sinking of the black hole or core stalling (Read et al. 2006; Petts et al. 2016).  \nWeconclude that the method of Tremmel et al. (2015) tends to sink the black holes too slowly, suggesting that the condition 𝑏 max = 𝜖 𝑔 is insufficient to account for unresolved dynamical friction. The method of Chen et al. (2022) tends to also underestimate the dynamical friction, due to the use of the Maxwellian approximation, particularly in the inner regions of an isotropic Hernquist halo and thus requires an enhanced value of 𝑏 max , unres to compensate for this underestimation. We confirm that 𝑏 max , unres = 6 𝜖 𝑔 is an effective correction at low resolution, but at high resolution it leads to a too fast sinking of black holes.', '4 A CALIBRATED DYNAMICAL FRICTION MODEL': "Aswehaveseenintheprevioussection,dynamicalfriction estimators like those of Tremmel et al. (2015) and Chen et al. (2022) are able to model the long-term orbital evolution of individual black holes, rather than sinking them immediately as in repositioning methods, and, with the right choice of parameters, these methods can sink black holes to the centres of their galaxies on a correct timescale. We will therefore adopt a similar form for our subgrid dynamical friction model:  \na df , M BH = -4 𝜋𝐺 2 𝜌 ( 𝑣 < 𝑣 BH ) 𝑣 3 𝑀 BH ln 𝜂𝜖 𝑔 𝑏 min v , (8)  \nwhere 𝜂 is a factor by which the softening is multiplied to obtain the maximum unresolved impact parameter. This is the value that we will attempt to calibrate through comparison with KETJU.  \nAnother choice we must make is the definition of the minimum impact parameter, 𝑏 min . Given the success of the Chandrasekhar analytical prediction in Section 2.2 compared to that of KETJU, we choose to define the minimum impact parameter as:  \n𝑏 min = 𝐺𝑀 BH ( 𝑣 2 BH + 2 3 𝜎 2 3D ) , (9)  \nwhich additionally sets the 'true' maximum impact parameter 𝑏 max = 𝑅 BH .  \nFigure 4. The contribution of the second-order diffusion coefficients to the change in the kinetic energy of the black hole, 1 2 𝐷 [ ( Δ 𝑣 ∥ ) 2 ] + 1 2 𝐷 [ ( Δ 𝑣 ⊥) 2 ] , compared to the contribution of dynamical friction, 𝑣𝐷 [ Δ 𝑣 ∥ ] . The ratio of energy contributions is computed as a function of time for the black hole orbits shown in the left panel of Fig. 1, but assuming unsoftened interactions. The terms are computed through evaluation of Eqns. (10), (11) and (12) for a Hernquist halo. The black dashed line shows where the contributions of the two terms are equal. The progressively darker shaded regions indicate where the contribution of stochasticity is above the 1 and 10 and 100 percent level of the dynamical friction contribution. Orbits of black holes with mass ratios as high as 𝑀 BH / 𝑀 DM = 100 can be substantially altered by stochastic effects (see also Fig. 11).  \n<!-- image -->", '4.1 Dynamical friction vs. stochasticity in collisional systems': 'Since we are attempting to reproduce the sinking of a black hole as it would occur in a smooth dark matter potential, where the masses of background particles tend to zero, we must be wary of how discreteness can affect black hole dynamics. Dynamical friction is one of the effects arising from gravitational encounters, which causes a black hole to lose its kinetic energy; however, scattering by massive background particles can also increase the kinetic energy of a black hole.  \nIn a collisionless system, the phase-space density, 𝑓 , around a test particle does not change with time. However, if gravitational collisions are present, the phase space density changes at a rate described bythe encounter operator Γ [ 𝑓 ] . One way to obtain the Chandrasekhar formula is via the diffusion coefficients of the Fokker-Planck approximation to the encounter operator (Spitzer 1987; Binney & Tremaine 1987). Namely, the first diffusion coefficient  \n𝐷 [ Δ 𝑣 ∥ ] = -16 𝜋 2 𝐺 2 ( 𝑀 BH + 𝑚 ∗ ) ln Λ 𝑣 2 ∫ 𝑣 0 d 𝑣 ∗ 𝑣 2 ∗ 𝑓 ( x , 𝑣 ∗ ) (10)  \nis just the dynamical friction formula of Eq. (1). The method of, e.g., Tremmel et al. 2015 neglects the additional \'cooling\' term arising from the mass 𝑚 ∗ of the background particles and thus assumes 𝑚 ∗ ≪ 𝑀 BH . The method of Damiano et al. (2024) includes the 𝑚 ∗ term in their correction for unresolved dynamical friction; however, when the black hole mass is lower than the dark matter particle mass,  \nthe resultant dynamical friction force is dominated by the 𝑚 ∗ term, which would otherwise be negligible in a well-resolved system 3 .  \nFor sufficiently large mass of the background particles compared to the black hole, the two second-order diffusion coefficients become important. The first second-order diffusion coefficient,  \n𝐷 [( Δ 𝑣 ∥ ) 2 ] = 32 𝜋 2 𝐺 2 𝑚 ∗ ln Λ 3 × " ∫ 𝑣 0 d 𝑣 ∗ 𝑣 4 ∗ 𝑣 3 𝑓 ( x , 𝑣 ∗ ) + ∫ ∞ 𝑣 d 𝑣 ∗ 𝑣 ∗ 𝑓 ( x , 𝑣 ∗ ) # , (11)  \ndescribes the stochastic contribution to the change of velocity squared per unit time in the direction of motion of the black hole. The second second-order coefficient,  \n𝐷 [( Δ 𝑣 ⊥ ) 2 ] = 32 𝜋 2 𝐺 2 𝑚 ∗ ln Λ 3 × " ∫ 𝑣 0 d 𝑣 ∗ 3 𝑣 2 ∗ 𝑣 -𝑣 4 ∗ 𝑣 3 ! 𝑓 ( x , 𝑣 ∗ ) + 2 ∫ ∞ 𝑣 d 𝑣 ∗ 𝑣 ∗ 𝑓 ( x , 𝑣 ∗ ) # , (12)  \ndescribes the heating in the direction perpendicular to the black hole motion 4 .  \nWe note that the order of magnitude of the second-order coefficients is determined primarily by the mass of the background particles, while the first-order coefficient (the strength of dynamical friction) depends on the black hole mass. Thus, the overall effect of stochastic heating on the sinking timescale of the black hole depends sensitively on the ratio 𝑀 BH / 𝑀 DM . The coefficient in Eq. (11) acts opposite to the direction of the dynamical friction force and accounts for the velocity dispersion of the slow-moving background particles and the typical velocity of the fast-moving particles.  \nIn Fig. 4 we show the contribution of the second-order diffusion coefficients to the energy change of the black hole compared to the contribution of the first-order (dynamical friction) coefficient. Overall, the stochastic terms become important for low 𝑀 BH / 𝑀 DM and, particularly, in regions of the halo where the velocity dispersion of background particles is larger than the black hole velocity. We note that even at 𝑀 BH / 𝑀 DM = 100, stochasticity can increase the kinetic energy of the black hole by a few per cent (at each timestep). Over time, these changes can add up and substantially alter the orbit of the black hole. It is therefore important that we are aware of the effects of the noise in the context of our model.  \nTherefore, in order to recover the dynamical friction force on a black hole in a smooth potential ( 𝑀 DM -→ 0) from the simulations, the full correction should take the form:  \na corr = a df , M BH -a df , m ∗ -a noise , ∥ -a noise , ⊥ , (13)  \nwith a noise , ∥ described by:  \na noise , ∥ = 4 𝜋𝐺 2 𝜌 ( 𝑣 < 𝑣 BH ) 𝑣 3 𝑚 ∗ 𝜎 2 3D 3 𝑣 2 ln 𝑅 𝜂𝜖 𝑔 v , (14)  \nwhere we have adopted the Maxwellian approximation in Eq. (11) (Binney & Tremaine 1987; Merritt 2013).  \nThrough explicit evaluation of Eqns. (10), (11) and (12), we have  \ndetermined that for an isotropic Hernquist halo the second-order diffusion term acting parallel to the direction of motion of the black hole, a noise , ∥ , approximately cancels out the additional cooling term, a df , m ∗ , and only dominates near the centre of the halo. The secondorder perpendicular term, a noise , ⊥ , is approximately equal in magnitude to the parallel term, suggesting that Eq. (14) is a good approximation to the magnitude of the noise terms (see Appendix A). However, we also found that the Maxwellian approximation typically overestimates the stochastic term in a Hernquist halo by a factor of ∼ 3, while the numerical evaluation of local 𝜎 3D / 𝑣 is extremely noisy and can lead to order-of-magnitude changes in the black hole acceleration. Moreover, it is unclear how one would account for the second-order terms perpendicular to the black hole motion, as we cannot at any given point determine the net direction of the stochastic force.  \nFor these reasons, we do not include the correction in Eq. (14) in this work. Instead, in the following we will derive an optimal value of 𝜂 that in effect marginalizes over these additional effects. For 𝑀 BH / 𝑀 DM > 100, the stochasticity contribution is at sub-per-cent level and will thus not influence our results (i.e. the optimal 𝑏 max , unres should only depend on the spatial resolution). Nevertheless, as we will discuss in Section 4.6, the stochastic effects can be substantial for 𝑀 BH / 𝑀 DM ≤ 100. We will therefore adopt an order-of-magnitude approximation for the overall effect of the second-order diffusion coefficients.', '4.2 Calibrating the maximum unresolved impact parameter': 'An obvious approach to calibrating 𝜂 would be to directly probe a range of values of this parameter by running the set of simulations from Table 1 with the corresponding dynamical friction correction. However, in this work we take a different approach that allows us to probe a wide range of 𝜂 at low computational cost. Our procedure is as follows:  \n(i) For each simulation with no dynamical friction prescription (left panel of Fig. 1), we take the position, velocity and acceleration of the black hole in the direction of motion, 𝑎 GADGET4 , at each timestep (excluding the times after black hole crosses within 2 kpc of the halo centre, where the density profile is suppressed compared to the Hernquist formula and where, at higher resolution, we would expect the local density centre of the halo to move in response to the motion of the black hole).  \n(ii) For a black hole inside an analytic Hernquist potential, with the same position and velocity vectors as the simulations in step (i), we compute the expected acceleration in the direction of motion of the black hole due to dynamical friction together with the effect of the external potential, 𝑎 analytical , following Eq. (1).  \n(iii) In a dark-matter-only simulation of the same halo and with the same gravitational softening, we measure the expected noise in the acceleration at the position of the black hole, 𝜎 GADGET4 , over multiple snapshots to account for noise due to the initial conditions. The noise arises here due to the discreteness of the 𝑁 -body potential and we assume it to be isotropic at any given location.  \n(iv) We sample values of 𝑏 max , unres in the range 0-100 kpc, and apply dynamical friction via Eq. (8), 𝑎 Eq . 13 . We assume that the density of particles moving slower than the black hole as well as the 3D velocity dispersion are known exactly. Using the affine invariant Markov Chain Monte Carlo ensemble sampler EMCEE (ForemanMackey et al. 2013), for each 𝑏 max , unres sample we compute the  \nFigure 5. EMCEE-generated samples of the distribution of optimal values of 𝜂 = 𝑏 max / 𝜖 𝑔 for the set of simulations in Table 1. The results are shown for the range 𝜂 = 0 -3 for clarity. No convergence is reached for 𝑀 BH / 𝑀 DM ≤ 10, though the peak of the distribution lies below 𝜂 ≈ 1. The optimal values of 𝜂 differ as the mass resolution and softening changes. Higher-resolution simulations prefer 𝜂 = 0 . 5 -0 . 8.  \n<!-- image -->  \nlog-likelihood:  \nln 𝐿 = -1 2 𝑡 = end ∑︁ 𝑡 = 0 ( 𝑎 GADGET4 + 𝑎 Eq . 13 -𝑎 analytical ) 2 𝜎 2 GADGET4 . (15)  \n(v) After a Monte Carlo Markov Chain (MCMC) burn-in length of 2500 steps with 250 chains (EMCEE walkers), we use the next 2500 steps to estimate the distribution of optimal 𝑏 max , unres choices.  \nWe note that, as an alternative method, we could instead determine the minimum resolved impact parameter 𝑏 min , res in simulations, given 𝑏 max = 𝑅 . This method is equivalent to our approach and we checked that this produces identical results. We additionally emphasise that the crucial aspect of this procedure is having a welldefined analytical 𝑏 min and 𝑏 max in the Coulomb logarithm. The former is well defined for a point particle as the impact parameter for a 90-degree deflection, while the latter is constrained through comparison with KETJU.  \nThe outcome of this exercise can be seen in Fig. 5, where we show the resulting distributions of 𝜂 = 𝑏 max , unres / 𝜖 𝑔 . It is clear from Fig. 5 that for the mass ratio 𝑀 BH / 𝑀 DM ≤ 10 our method has failed to converge on an optimal 𝜂 , although already at 𝑀 BH / 𝑀 DM = 5 we can see that this value needs to be above zero. This result is not surprising, as at the lowest resolution, 𝑀 BH / 𝑀 DM = 1, we do not expect to see any effects of dynamical friction. In this case, the orbit of the black hole is quite noisy, and in some instances the dynamical friction correction would need to be either small or large, and positive or negative, to obtain the acceleration that corresponds to the analytical prediction. This is still partially the case until 𝑀 BH / 𝑀 DM > 10. At 𝑀 BH / 𝑀 DM ≥ 100, however, we are already able to effectively narrow down the range of optimal values of 𝜂 , which appears to be within 0 . 5 -0 . 7.  \nWe also note that in our approach of examining the simulations from the left panel of Fig. 1, some of the black holes (particularly at the three lowest resolution levels) never approach the centre of the halo. This could suggest that the lack of a constraint on 𝜂 at this resolution arises due to this missing data; however, we have explored the radial dependence of our constraint (i.e. repeating the same procedure as outlined above but for radial bins of 2 kpc at  \nall resolution levels) and found consistent best values of 𝜂 down to the inner ∼ 2 𝜖 𝑔 , where the density profiles of 𝑁 -body haloes would naturally deviate from analytical predictions.', '4.3 Probing different values of the softening': "Our experiment in Fig. 5 was limited to the range of 'optimal' softening choices according to Power et al. (2003), and it is not yet clear whether the required optimal 𝜂 would be the same if the softening was chosen differently. We therefore repeat the procedure outlined in the previous section, now varying the softening for each background particle mass resolution. The softening values we use at each resolution level are the same as those in Table 1.  \nThe result of this is given in Fig. 6, where on the left we show the maximumunresolvedimpactparameterasafunctionofthesoftening, and on the right we show the optimal 𝜂 obtained through our method in each case. Firstly, we see that for all gravitational softening values there is no constraint on 𝑏 max , unres for the lowest resolution run, with 𝑀 BH = 𝑀 DM . At high resolution ( 𝑀 BH / 𝑀 DM ≥ 50 and 𝜖 𝑔 ≤ 2) there appears to be a clear correlation between the softening length and unresolved 𝑏 max , as well as a convergence around the optimal value of 𝜂 . Above 𝜖 𝑔 = 2 kpc, the relation appears to steepen, however within the error bars it is still consistent with the linear relation we see for higher-resolution simulations. The steepening of the relation is likely reflective of the unresolved central densities of the simulated halos, thus requiring stronger dynamical friction to sink the black holes consistently with the analytical prediction.  \nThe results for 𝑀 BH / 𝑀 DM = 5 and 𝑀 BH / 𝑀 DM = 10 have a particularly interesting shape, where at larger values of softening the relation is followed, while at values of the softening that are too low for this resolution, 𝜖 𝑔 < 2 kpc, the maximum unresolved 𝑏 max appears to increase. We interpret this as the result of enhanced two-body interactions that introduce additional noise in the black hole orbit. Nevertheless, within the errors, even these results are consistent with the general relation seen at higher resolution. We thus see no evidence of a relation between the mass resolution and the maximum unresolved values of 𝑏 max , which must be set only by the gravitational softening.  \nTo obtain the optimal value of 𝜂 , we fit a quadratic relation to the left panel of Fig. 7. This is shown with a solid grey line and corresponds to the best fit value  \n𝜂 optimal = min [ 0 . 46 ( 𝜖 g / kpc ) 0 . 17 10 0 . 11 ln 2 ( 𝜖 g / kpc ) , 1 . 0 ] , (16)  \nwhere we fixed the maximum value of 𝜂 to 1. The relation 𝑏 max = 𝜖 𝑔 is shown with the grey dashed line. We note that the definition of 𝑏 min of Tremmel et al. (2015) and Chen et al. (2022) is different to ours, so our relation would not hold in that case; however, with 𝑏 min = 𝐺𝑀 / 𝑣 2 BH , and through an approximation 𝑣 2 BH ≈ 1 3 𝜎 2 3D , we can deduce that an optimal value is 𝜂 ≈ 2 for Tremmel et al. (2015) and Chen et al. (2022). We show the optimal value of 𝜂 with this definition of 𝑏 min in Appendix B.  \nFinally, we wish to verify that our relation is not specific to the particular black hole mass and its corresponding orbit within our mock halo. To do this, we now repeat our test, while varying the halo and black hole properties. The result of this can be seen in Fig. 7, where we show constraints for the optimal vale of 𝜂 for a 10 8 M ⊙ black hole in a 10 12 M ⊙ halo, and a 10 7 M ⊙ black hole in a 10 11 M ⊙ halo. It can be seen that the obtained values agree well with our previously derived relation when the softening is normalized by the Hernquist scale radius 𝑎 HQ , suggesting that the variation in optimal 𝜂 is dependent on the rate of encounters experienced by the black hole in the simulation.  \n<!-- image -->  \nFigure 6. Left: the maximum unresolved impact parameter, 𝑏 max , unres , as a function of gravitational softening, 𝜖 𝑔 , obtained through the method described in Section 4.2. Different colours correspond to the mass resolution of a 10 13 M ⊙ halo hosting a 10 9 M ⊙ black hole and the respective black hole to dark matter particle mass ratio. For each of the 9 levels of mass resolution, simulations of the 10 9 M ⊙ black hole starting off at 20 kpc on a circular orbit were run for 9 different values of the softening. The circles indicate the median of the posterior distribution of the maximum unresolved impact parameter obtained with EMCEE for each simulation, and error bars show the 16 th and 84 th percentiles. The solid grey curve shows the least-squares fit of a quadratic relation to all data points with M BH / 𝑀 DM > 1. The dashed grey line shows 𝑏 max , unres = 𝜖 𝑔 (Tremmel et al. 2015). The arrows indicate the optimal softening at each mass resolution, according to Power et al. (2003). Right: same as the left panel, but for the constrained optimal 𝜂 = 𝑏 max , unres / 𝜖 𝑔 .  \n<!-- image -->  \nFigure 7. Similar to Fig. 6, but showing our results for black holes and haloes of varying mass. The lines show the optimal 𝜂 at each mass resolution for our fiducial 10 13 M ⊙ halo as a function of softening, normalized by the scale length of the Hernquist halo ( 𝑎 HQ = 80 kpc). The colours denote the mass resolution and arrows give the optimal softening at a given mass resolution for the 10 13 M ⊙ halo. Square and star symbols show derived values of 𝜂 at optimal softening for a 10 12 M ⊙ halo (10 8 M ⊙ black hole starting at 10 kpc) and a 10 11 M ⊙ halo (10 7 M ⊙ black hole starting at 5 kpc), respectively.  \n<!-- image -->  \nThe generalized relation for any given Hernquist-like halo may therefore be expressed as:  \n𝜂 optimal = min [ 103 . 7 ( 𝜖 g / a HQ ) 1 . 3 10 0 . 11 ln 2 ( 𝜖 g / a HQ ) , 1 . 0 ] . (17)  \nWe note that while 𝑎 HQ can be estimated with an on-the-fly halo finding algorithm, in practice, many cosmological simulations will not reach the softening of 10 -3 𝑎 HQ and will thus fall on the 'flat' part of our relations with 𝜂 ≈ 0 . 45 𝜖 𝑔 .", '4.4 The meaning of the optimal 𝜂': 'It may at first seem surprising that our derived optimal values of 𝜂 are smaller than 2.8 𝜖 𝑔 , given that this is the radius where the forces become fully Newtonian in GADGET-4 (i.e. particle scattering follows the full Newtonian case). To gain an intuitive understanding of this apparent discrepancy, we consider 2-body gravitational scattering. The use of the gravitational softening introduces a maximum scattering angle in a given interaction. The impact parameters that correspond to stronger particle deflections are indeed unresolved (for the spline softening we infer that the maximum deflection occurs at 𝑏 ≃ 1 . 4 𝜖 𝑔 , see White 1976; Ludlow et al. 2019b). However, this also implies that there is now an excess of particles with weaker deflections, effectively increasing the rate of encounters with impact parameters 𝑏 > 1 . 4 𝜖 𝑔 (though the exact dynamical friction effect will depend on the momentum exchange resulting from the scattering). Therefore, the correction to dynamical friction that one needs to apply should not only compensate for the unresolved impact parameters but also take away the excess encounters with weaker deflections, which is what our relation achieves.  \nA seemingly puzzling feature of our recovered relation is that, while the maximum unresolved impact parameter 𝑏 max , unres does decrease with the softening, it appears to decrease slower towards  \nsmaller values of the softening, resulting in a 𝜂 ( 𝜖 𝑔 ) relation that has a parabolic shape. However, we can also obtain a similar parabolic relation using equation 6 of White (1976), which computes the minimum impact parameter for an interaction of an extended object (e.g. a globular cluster) with a field star, where for the mass distribution of the subject we assume the distribution implied by the spline softening in GADGET-4. We note however, that the treatment of White (1976) assumes undeflected trajectories of field stars when computing their velocity change due to the gravitational encounter with an extended object and thus only sets an upper limit on their minimum impact parameter ( 𝑏 min ≤ 1 . 4 𝜖 𝑔 ).', '4.5 Optimizing the kernel size': 'As pointed out previously, our calibration procedure is based on the assumption that the density 𝜌 ( 𝑣 < 𝑣 BH ) is known exactly and so is the velocity dispersion, 𝜎 3D (which determines 𝑏 min ). However, as we have seen in Fig. 3, the density estimate is likely to be noisy. In our method, we should thus aim to minimize this noise while simultaneously maintaining the accuracy of the density and velocity dispersion estimates.  \nIn order to obtain an optimal kernel size at each resolution level, we first take a fixed black hole orbit from the highest-resolution run shown in the left panel of Fig. 2 (10 7 particles, 𝑀 BH / 𝑀 DM = 1000), defined by 3D positions and velocities, ( r BH ( 𝑡 ) , v BH ( 𝑡 ) ). We then compare the values of 𝜌 ( 𝑣 < 𝑣 BH ) and 𝜎 3D measured in dark-matteronly simulations for each ( r BH ( 𝑡 ) , v BH ( 𝑡 ) ) at each resolution level to those predicted by the Hernquist distribution function 5 . Note, however, that we opt to compare the cumulative density encountered by the black hole along its trajectory, " 𝜌 = " 𝑡 = current 𝑡 = 0 𝜌 ( 𝑣 < 𝑣 BH ) , with the analytical prediction. This is because a density estimate that is noisy but on average accurate may still sink the black hole on the correct timescale, compared to a smooth estimate with a smaller sum of individual density residuals, but a position-dependent bias. We additionally estimate the errors for each " 𝜌 and 𝜎 3D by incrementally rotating the trajectory of the black hole in 3 dimensions, such that we have a distribution of " 𝜌 and 𝜎 3D at fixed ( 𝑑 BH ( 𝑡 ) , v BH ( 𝑡 ) ), where 𝑑 BH is the distance of the black hole from the centre of the halo.  \nWe note that as the black hole approaches the halo centre and encounters steeper density gradients, a large kernel size can bias the density estimate. For this reason, we define a \'maximum resolved density\' 6 , 𝜌 max ≈ 𝜌 200 , crit √︁ 𝑁 p . We found that this density approximates well the transition in the density profile where smaller kernel sizes are more advantageous than larger ones. We then sample the kernel sizes in the range 10 to 1000 particles using EMCEE for the trajectory of the black hole in regions of the halo where the local density is smaller than 𝜌 max , computing the log-likelihood:  \nln 𝐿 𝑘 = -1 2 𝑡 = 𝑡 ( 𝜌 max ) ∑︁ 𝑡 = 0 GLYPH<18> " 𝜌 ( 𝑁 𝑘 ) -" 𝜌 analytical 𝜎 " 𝜌 ( 𝑁 𝑘 ) + 𝜎 3D ( 𝑁 𝑘 ) -𝜎 3D , analytical 𝜎 𝜎 3D ( 𝑁 𝑘 ) GLYPH<19> , (18)  \nwhere 𝜎 " 𝜌 and 𝜎 𝜎 3 𝐷 are the errors on the estimates.  \nThe results can be seen on the left panel of Fig. 8, through which we derive a relation:  \nlog 10 𝑁 k =          min ( 0 . 36 log 10 𝑁 p + 0 . 33 , 3 ) if 𝜌 < 𝜌 max , log 10 64 if 𝜌 > 𝜌 max , log 10 64 if 𝜖 𝑔 / 𝑎 HQ > 0 . 03 , (19)  \nwhere 𝑁 p is the number of particles in the halo, and 𝑁 k is the optimal kernel size. We chose to limit the number of particles in the kernel to a minimum of 64 and a maximum of 1000. The middle panel of Fig. 8 shows how the accuracy of estimates of " 𝜌 ( 𝑣 < 𝑣 BH ) using our optimal kernel size compare to the values predicted by the Hernquist distribution function. The thin lines shows the result for a fixed kernel size of 64 particles as in Tremmel et al. (2015), where the benefit of our approach at higher resolution is particularly evident. In practice, as a typical treatment for black hole seeding in cosmological simulations involves identifying the mass of the host halo, 𝑁 p can be approximated simply as 𝑀 halo / 𝑚 DM .', '4.6 A treatment for low-mass black holes': 'Finally, we would like our method to be able to deal with the case of low-mass black holes (i.e. black holes with mass less or equal to the dark matter particle mass). As we have already seen, for the case 𝑀 BH / 𝑀 DM = 1 the black hole behaves in a similar way as the dark matter and does not sink. For 𝑀 BH / 𝑀 DM < 1, we furthermore expect to see effects of mass segregation (Ludlow et al. 2019a), where lower-mass particles are scattered to larger radii by higher-mass ones.  \nWithin the framework of our model, we can make low-mass black holes sink either by increasing their dynamical mass, increasing the magnitude of the Coulomb logarithm, or fully accounting for the heating due to stochasticity for a given background particle mass.  \nWhile increasing the dynamical mass of the black hole may be in part physically motivated as many black holes are observed to live inside star clusters (see Kormendy & Ho 2013; Neumayer et al. 2020 and references therein), we do not at present have a subgrid prescription for the formation of these clusters. It is also unclear to what extent would these clusters enhance the dynamical mass of the black hole, with observations suggesting star clusters both more and less massive than the black holes they host, and whether this boost in mass would guarantee sinking at given mass resolution.  \nIn Fig 9, in black, we show the circular orbit of a 10 8 M ⊙ and a 10 7 M ⊙ black hole in a 10 13 M ⊙ halo. In both cases, 𝑀 BH / 𝑀 DM = 1. Coloured lines show the orbit of a black hole as we increase its dynamical mass by a factor 𝑘 between 1 and 10. It can be seen that the mass ratios 𝑀 dyn / 𝑀 DM < 4 are often not enough to guarantee sinking of the black hole without added dynamical friction (Bellovary et al. 2021). On the other hand, 𝑀 dyn / 𝑀 DM ≳ 4 are often sufficient to sink the black hole. However, while increasing the dynamical mass may be physically justified if the black hole is embedded in a more massive star cluster, it may also lead to undesirable effects such as artificial scattering of background particles by the black hole (Wurster & Thacker 2013). In any case, Fig. 9 demonstrates that there is no particular dynamical mass scaling that would reproduce the black hole orbit for all initial conditions.  \nOur calibration method, described in Section 4.2 with results shown in Figs. 6 and 7, suggests a strategy for increasing the magnitude of the Coulomb logarithm to compensate for the stochasticity, as our derived optimal values of 𝜂 for each 𝑀 BH / 𝑀 DM marginalise over the impact of noise. However, 𝜂 is not well constrained for low mass ratios and is close to a flat distribution, as shown in Fig. 5, due to  \n<!-- image -->  \n<!-- image -->  \nFigure 8. Left: the number of particles in the kernel required to obtain an accurate estimate of " 𝜌 ( 𝑣 < 𝑣 BH ) as a function of the number of particles in the halo. The points are coloured by the size of the gravitational softening from largest (blue) to highest (orange). The error bars are shown in each case, but are in many cases smaller than the point size. The black line is the best linear fit to the relation. The maximum number of particles in the kernel is limited to 1000 (grey dashed line) and the minimum is capped at 64 particles (grey shaded region). Middle: histograms with solid lines show the accuracy of the cumulative density encountered by the black hole in halos with varying mass resolution with a fixed softening of 1.4 kpc. The densities were computed using our optimal-kernel-size method and compared to the prediction from the Hernquist distribution function. Thin dashed lines show the estimate using a kernel defined by 64 nearest neighbouring particles as in Tremmel et al. (2015). Right: the accuracy of the estimate of the 3D velocity dispersion with a fixed softening of 1.4 kpc. Solid lines show the results using our optimal kernel size, and dashed lines show the results for a kernel with 64 particles.  \n<!-- image -->  \n<!-- image -->  \nFigure 9. Left: the effect of increasing the dynamical mass of black holes on their orbital trajectories. In black we show the analytical prediction for the orbit of a 10 8 M ⊙ black hole in a 10 13 M ⊙ halo, starting on a circular orbit at 10 kpc. The initial orbit with 𝑀 BH / 𝑀 DM = 1 is shown with a thin dark blue line. The dynamical mass of the black hole is then scaled by a factor 𝑘 in the range 1 to 10 times the dark matter particle mass, resulting in faster black hole sinking. No dynamical friction correction is applied. The thick red line shows the result of setting the dynamical mass of the black hole to 𝑀 dyn = 𝑀 BH + 𝑀 DM and applying the correction for unresolved dynamical friction via Eqn. (8). Right: same as the left panel, but for a 10 7 M ⊙ black hole starting at 5 kpc. The black hole trajectories have been smoothed over 0.1 Gyr for clarity.  \n<!-- image -->  \nstochastic forces having comparable magnitude to that of dynamical friction i.e at any given position the force that the black hole experiences can deviate significantly from that due to a smooth Hernquist halo. We thus conclude it would be more advantageous to derive an estimator for the overall effect of stochasticity from first principles.  \nAs discussed in Section 4.1, while we can obtain an estimate of the magnitude of the stochastic acceleration, we cannot at any given point determine the net direction of the stochastic kicks. In  \nsome cases, stochastic effects can ultimately lead to a more rapid black hole sinking (right panel of Fig. 1). However, as the effect of stochasticity is to increase the kinetic energy of the black hole, we can derive an energy-loss term to counteract this.  \nWethus will assume that the overall effect of the second-order diffusion coefficients, working against the effect of dynamical friction,  \nhas the form:  \n𝑎 noise = 4 𝜋𝐺 2 𝜌 ( 𝑣 < 𝑣 BH ) 𝑀 DM 𝑣 3 ln GLYPH<18> 𝑅 𝜂𝜖 𝑔 GLYPH<19> v if 𝑅 > 𝜂𝜖 𝑔 , (20)  \nwhere 𝑀 DM is the background particle mass, which brings the noise correction to the right order of magnitude (see Appendix A). Since particle scattering below 𝜂𝜖 𝑔 is unresolved, only impact parameters above 𝜂𝜖 𝑔 contribute, giving the Coulomb logarithm ln Λ = 𝑅 𝜂𝜖 𝑔 . We note that this is different to applying the term a df , m ∗ as in the work of Damiano et al. (2024) - their additional \'cooling\' force accounts for unresolved dynamical friction from the mass of background particles (assuming their masses are the true masses of background stars and dark matter), while Eq. (20) attempts to negate the resolved heating of the black hole by massive background particles.  \nWenote that the application of the noise correction is, in principle, equivalent to setting the dynamical mass of the black hole in the gravity solver to  \n𝑀 dyn = 𝑀 BH + 𝑀 DM , (21)  \nwhile the correction for the unresolved part of dynamical friction is still computed for a black hole mass 𝑀 BH . The mass rescaling approach may also be more advantageous in practice, as computing the distance to the centre of the halo at each time step entails a higher computational cost. Additionally, the definition 𝑏 max = 𝑅 BH comes from the assumption of a Hernquist halo (Just et al. 2011) and can change for other density distributions. The exact definition of 𝑏 max is irrelevant in mass boosting. Nevertheless, in practice, we find that directly applying the noise term will sink the black holes more aggressively and guarantees sinking, while increasing the dynamical mass as per Eq. (21) reproduces the orbital path of the black holes more accurately on average, but in some cases can result in black hole sinking that is substantially delayed or even suppressed.  \nTheresult of the application of this stochasticity correction through a rescaling of the dynamical mass, together with the correction for unresolved dynamical friction, is shown with a thick red line in Fig 9. In neither case does the correction fully reproduce the true trajectory of the black hole, however, in both cases the black holes sink to the centre of the galaxy without the need for a substantial increase of their dynamical mass.', '4.7 Black holes in the regime 1 ≤ 𝑀 BH / 𝑀 DM ≤ 5': 'We now explore how our stochasticity correction behaves in the regime 1 ≤ 𝑀 BH / 𝑀 DM ≤ 5. This is shown in Fig. 10, where for a fixed halo and black hole mass, we gradually increase the mass resolution such that 𝑀 BH / 𝑀 DM varies from 1 to 5. The dark lines show the median trajectory of 6 black hole seeds for each mass ratio. We apply the stochasticity correction via Eq. (20). The correction clearly tends, on average, to sink the black holes too fast, with the effect being strongest for lower 𝑀 BH / 𝑀 DM . The error bars range from the earliest to the latest sinking times of the black holes, where the mass ratio appears to set the earliest black hole sinking time, while the maximum sinking time does not seem to depend strongly on the resolution. The light lines show individual realisations of black holes at each mass ratio where only the dynamical friction correction is applied. In some cases, this faithfully recovers the black hole orbit ( 𝑀 BH / 𝑀 DM = 2, 2 . 5, and 4 . 5), but for others the black holes do not sink within the age of the Universe. Interestingly, this does not seem to depend strongly on the mass ratio in this regime (some black holes with 𝑀 BH / 𝑀 DM = 4 and 5 do not sink even with a dynamical friction correction). We conclude that a stochasticity correction is required to ensure that a black hole sinks to the centre of a simulated halo,  \nFigure 10. The results of applying our dynamical friction correction, coupled with the correction for stochasticity, to a 10 7 M ⊙ black hole in a 10 13 M ⊙ halo starting at 4 kpc. Darker lines show median black hole trajectories for different 𝑀 BH / 𝑀 DM ratios as the halo mass resolution changes. The error bars show the range of sinking times of black holes from the sample of 6 initial conditions used to compute the median trajectories. The lighter lines show examples of trajectories where only a dynamical friction correction is applied, i.e. no stochasticity correction. The black hole trajectories have been smoothed over 0.15 Gyr for clarity.  \n<!-- image -->  \nhowever the inherent randomness of this process makes it difficult to recover black hole orbits accurately for any individual case in the regime 𝑀 BH / 𝑀 DM ≤ 5.  \nWe have repeated this test with rescaling the dynamical mass of the black hole to 𝑀 BH + 𝑀 DM instead of applying the stochasticity correction directly. We found that up to 𝑀 BH / 𝑀 DM ≤ 5, 10 per-cent of the initial conditions lead to substantailly delayed sinking of the black holes (the black holes do not sink before 13 Gyr) and 13 per cent of initial conditions show no or little sign of sinking over 13 Gyr, while the median sinking times for each 𝑀 BH / 𝑀 DM lie between 8 and 12 Gyr. Directly applying the stochastic correction via Eqn. (20), on the other hand, at all times guarantees the sinking of black holes.  \nWebelieve the discrepancy between the two arises for two reasons. One is the scattering that the black hole of higher dynamical mass induces, which can alter the properties of the surrounding medium. Another is that, when the correction is applied directly the net force is always exactly opposite to the direction of motion of the black hole, maximising the reduction in kinetic energy.  \nIn conclusion, in the regime 1 ≤ 𝑀 BH / 𝑀 DM ≤ 5, the predictions for individual black hole trajectories become highly uncertain when using a dynamical friction subgrid model such as the one described in this work. To guarantee black hole sinking, one needs to account for the effects of stochasticity e.g. via Eq. (20). The dynamical mass rescaling as in Eq. (21) ( 𝑀 dyn = 𝑀 BH + 𝑀 DM ) works well on average, but in ∼ 13 per cent of cases this mass scaling is too weak to sink the black hole. Increasing the dynamical mass to 𝑀 dyn ≥ 4 𝑀 DM would guarantee sinking, but at the price of lost accuracy in the orbital trajectory. We also note that the regime of 𝑀 BH / 𝑀 DM < 1 has not  \nbeen explored in this work. As such, in the framework of this model, we recommend that the black holes are seeded at least with the mass of the dark matter particles.', '5 APPLYING THE MODEL': 'In this Section, we briefly summarise the necessary components of our model and the requirements for its implementation in a cosmological simulation. We then apply the model to a number of test cases.', '5.1 Summary of the model': "- (i) When the black hole is seeded in a halo of given mass, the number of particles in the halo is approximated as 𝑁 𝑝 = 𝑀 halo / 𝑀 DM and the 'maximum resolved density' as 𝜌 max = 𝜌 200 , crit √︁ 𝑁 𝑝 . The kernel size for estimating the local density and velocity dispersion for each black hole can then be determined through Eq. (19). This value can be updated when the halo finder is applied again 7 .\n- (ii) The scale radius of each halo, 𝑎 HQ , may be estimated via, e.g., the mass-concentration relation in Λ CDM (Ludlow et al. 2016) and eq. 2 of Springel et al. (2005a) 8 . For a given gravitational softening of the black hole, 𝜖 𝑔 , Eq. (17) is used to determine the maximum unresolved impact parameter, 𝑏 max , unres = 𝜂 ( 𝜖 𝑔 ) 𝜖 𝑔 . We emphasise that the relation in Eq. (17) is only applicable to codes that use the same spline softening as in GADGET and AREPO (Springel 2010; Springel et al. 2021).\n- (iii) The black hole is assigned a dynamical mass 𝑀 dyn = 𝑀 BH + 𝑀 DM , which is the mass used in the gravity solver. The dynamical mass should be updated as long as M BH / 𝑀 DM ≤ 100, subject to the following constraints:\n- · The black hole dynamical mass should always be at least twice the dark matter particle mass 9 .\n- · As an alternative to dynamical mass rescaling, one may wish to apply the noise correction directly to the black hole, in which case one needs to estimate at each timestep the distance of the black hole from the centre of the halo 10 , 𝑅 .\n- (iv) At each timestep, the local density of particles moving slower than the black hole, 𝜌 ( 𝑣 < 𝑣 BH ) , and the 3D velocity dispersion, 𝜎 3D , are computed from the properties of the nearest 𝑁 𝑘 particles. The minimum impact parameter is computed as 𝑏 min = 𝐺𝑀 BH /( 𝑣 2 BH + 2 / 3 𝜎 2 3D ) , where the black hole velocity is relative to the mean velocity of the surrounding particles. In case of multiple populations (e.g. stars and dark matter), the properties are computed separately for each population and so is the resultant  \ndynamical friction. 𝑁 𝑘 should also be computed for different populations. If a negligible density of particles is present in a given population, the dynamical friction contribution may be ignored.  \n- (v) The correction for unresolved dynamical friction is applied using Eq. (8), with the Coulomb logarithm computed as ln Λ = 𝑏 max , unres / 𝑏 min . The correction is applied as long as 𝑏 min < 𝑏 max , unres . If the dynamical mass of the black hole has not been rescaled, a stochasticity correction should be applied using Eq. (20).", '5.2 Validation': 'We will now demonstrate how our method performs in a variety of setups. In particular, we will test our fiducial circular orbit setup, a more eccentric orbit and the test case presented in Tremmel et al. (2015).', '5.2.1 Fiducial set-up': 'In Fig. 11, we show the median black hole trajectories generated from 6 random realisations of the Hernquist halo at 7 different mass resolutions. The top panel shows simulations with the optimal Power et al. (2003) softening, while the bottom panel shows the results with a fixed softening of 0.44 kpc. The left and right panels give the results without and with the stochasticity correction, respectively. We compare our results to those using the Tremmel et al. (2015) and Chen et al. (2022) models for the highest-resolution case ( 𝑁 p = 10 7 , 𝑀 BH / 𝑀 DM = 1000).  \nIt is evident that our calibrated formula reproduces well the sinking time of the black holes from the analytical and KETJU predictions. With 𝑀 BH / 𝑀 DM ≥ 100 the sinking times are well converged, while at lower mass ratios a stochasticity correction is clearly required, as deviations from the prediction become apparent even at large radii. The importance of optimal softening choice is especially evident in the lower panel of Fig. 11, where the softening of 𝜖 𝑔 = 0 . 44 kpc is clearly too small for halos with 𝑁 p ≤ 10 5 M ⊙ , leading to strong stochastic effects. Black holes with 𝑀 BH / 𝑀 DM = 1, for instance, tend to get ejected away from the halo centre, even with application of the unresolved dynamical friction force.  \nThe fluctuations seen in the median black hole orbit at 11 Gyr for cases with 𝑀 BH / 𝑀 DM = 10 and 𝑀 BH / 𝑀 DM = 50 correspond to individual black holes that only arrive to the centre of the halo at this time. In one case with 𝑀 BH / 𝑀 DM = 10, the black hole sinks to the centre of the halo at 6 Gyr, but then gets ejected from the centre before sinking again (a rise in the median profile between 6 and 8 Gyr). The stochasticity correction (Eq. 20) avoids these effects (bottom right panel). For the case with 10 4 particles, a constant density core on the scale of ∼ 7 kpc develops, leading to the black hole stalling around this radius and additionally making it difficult to define the halo density centre.', '5.2.2 An eccentric orbit': 'We now test our method on a more eccentric black hole orbit. The 10 9 M ⊙ black hole in a 10 13 M ⊙ halo has a starting velocity 𝑣 BH = 0 . 5 𝑣 circ at 20 kpc, resulting in a sinking time of 2 . 3 -2 . 4 Gyr. The results of the test are shown in Fig. 12. This time, we see a mild difference in sinking time ( ∼ 0 . 15 Gyr) between the analytical prediction and KETJU. This likely stems from the limited resolution at the centre of the halo. The thick lines show the median black hole trajectories, where we see that the dynamical friction and stochastic  \nFigure 11. Top: the median trajectories of 10 9 M ⊙ black holes in a 10 13 M ⊙ halo computed using our calibrated dynamical friction model. The results are shown for 7 mass resolutions and their corresponding optimal softening values (see Table 1). The grey dotted and dashed lines show the sinking time of the black hole at the highest background particle mass resolution using the Chen et al. (2022) and Tremmel et al. (2015) models, respectively. The error bars show the range of sinking times computed from 6 simulations at each resolution with differing random seeds in the initial conditions. The sinking time is defined as the time when the black hole first comes one softening length away from the centre of the halo. The left and right panels show the results with and without the application of the noise correction, respectively. Bottom: same as the top panel, but using a fixed softening of 𝜖 𝑔 = 0 . 44 kpc.  \n<!-- image -->  \ncorrection are in this case slightly too weak, with black hole sinking time delayed on average by ∼ 0 . 2 Gyr. This may seem contradictory to our findings in Figs. 10 and 11, but we emphasize that both of those represent extreme cases with long black hole sinking timescales over which stochastic effects become obvious.', '5.2.3 An example from Tremmel et al. (2015)': 'Finally, we apply our model to the initial conditions similar to the H/i.pc/g.pc/h.pcR/e.pc/s.pc test case presented in Tremmel et al. (2015). This in-  \ncludes a 𝑀 200 = 2 × 10 11 M ⊙ halo with concentration 𝑐 = 4 . 5 and 𝑅 200 = 115 kpc, in which a black hole with mass 10 6 M ⊙ starts off at 2 kpc, with a velocity 𝑣 = 0 . 1 𝑣 circ . The dark matter particle mass in our simulations is 2 . 7 × 10 3 M ⊙ , resulting in 𝑀 BH / 𝑀 DM ≈ 380. The gravitational softening is set to 𝜖 𝑔 = 77 pc, as in Tremmel et al. (2015). In Fig. 13, the red line shows the black hole trajectory using our dynamical friction model, and the blue line shows the trajectory of the black hole using the Tremmel et al. (2015) model. The sinking time prediction for the analytical model and KETJU differs by ∼ 0.05 Gyr. This could be the result of minor stochastic effects we expect for mass ratios 100 < 𝑀 BH / 𝑀 DM < 500 (see Fig. 4).  \nFigure 12. The median trajectories of 10 9 M ⊙ black holes set on an eccentric orbit in a 10 13 M ⊙ halo. The median trajectories are computed from 6 random seeds at each resolution level (same as Fig. 11). The error bars show the range between the earliest and the latest black hole sinking time. The black hole dynamics were computed using our correction for the unresolved dynamical friction and the correction for stochasticity. The grey dotted and dashed lines show the sinking time for 𝑀 BH / 𝑀 DM = 1000 predicted by the Chen et al. (2022) and Tremmel et al. (2015) models, respectively.  \n<!-- image -->  \nFigure 13. The trajectory of a 10 6 M ⊙ black hole in a 2 × 10 11 M ⊙ halo from the H/i.pc/g.pc/h.pcR/e.pc/s.pc example presented in Tremmel et al. (2015). The red line shows the result of applying our model, while the blue line is the result of the Tremmel et al. (2015) model. The black solid and dashed lines are the analytical and KETJU predictions, respectively. All simulations were carried out with 𝑀 BH / 𝑀 DM = 376.  \n<!-- image -->  \nWe also note that this was the most computationally expensive of our simulations in this work, as the accuracy of the time integration in the analytical model (and consequently the simulations with our model and the model of Tremmel et al. 2015) had to be increased to 0.1 per cent in order to accurately capture the eccentric orbit of the black hole, as confirmed by convergence tests. Our model follows the analytical prediction and that of KETJU almost exactly, while the model of Tremmel et al. (2015) results in the black hole sinking time delayed by almost 0 . 2 Gyr.', '6 CONCLUSIONS': "Dynamical friction is one of the main factors governing black hole motion, driving black holes to the centres of galaxies where they can accrete and grow. It is therefore essential to model this process accurately in cosmological hydrodynamics simulations. Unfortunately, the spatial and mass resolution required to capture the full extent of dynamical friction comes at an extremely high computational cost, especially for large cosmological volumes. In particular, the use of the gravitational softening limits the scattering of local particles by the black hole, which would otherwise contribute a significant fraction of the total dynamical friction force. Coarse-grained representation of background stars and dark matter also introduces artificial heating of the black hole orbit. The latter effect is particularly problematic in cases where the black hole particle mass is of similar magnitude to that of the background stars and dark matter.  \nA number of numerical models aimed at compensating for the unresolved part of the dynamical friction force have been introduced in the literature (Tremmel et al. 2015; Pfister et al. 2019; Chen et al. 2022; Damiano et al. 2024). These models assume that the effect of dynamical friction in an 𝑁 -body simulation is well-described by the Chandrasekhar formula (Chandrasekhar 1943) in which the minimum impact parameter of the Coulomb logarithm, ln Λ = 𝑏 max 𝑏 min , is determined by the spatial resolution of the simulation ( 𝑏 min ≈ 𝜖 𝑔 ). They then apply a correction for the unresolved part of dynamical friction, with ln Λ unresolved = 𝑏 max , unres 𝑏 90 ≈ 𝜖 𝑔 𝑏 90 , where 𝑏 90 is the impact parameter for a 90-degree deflection. Such a correction has been shown to effectively sink black holes to galactic centres on roughly correct timescales. In the case of low-mass black holes ( 𝑀 BH ∼ 𝑀 DM ), these models artificially boost the dynamical mass of the black hole or increase the strength of dynamical friction to overpower the effects of the numerical noise, with the main goal of keeping low-mass black holes near the centres of their host galaxies.  \nNevertheless, it has remained unclear what is the exact spatial scale up to which dynamical friction is unresolved in 𝑁 -body simulations and how this scale relates to the choice of gravitational softening. The effect of numerical noise on black hole motion has also been largely ignored in the literature. These are the issues we have aimed to address in this work. By comparing the dynamics of black holes in simulations with varying spatial and mass resolution to high-resolution simulations run with the code KETJU, in which the gravitational interaction of black holes are not softened, we were able to calibrate a subgrid model based on the Chandrasekhar formula to accurately follow orbits of black holes and their sinking to halo centres. In particular:  \n(i) We have shown that the definition of the impact parameter for a 90-degree deflection, 𝑏 min , is important. The definition 𝑏 min = 𝐺𝑀 BH /( 𝑣 2 BH + 2 / 3 𝜎 2 3 𝐷 ) , taking into account random motions in the typical velocity of an encounter, ensures that 𝑏 min is generally converged when changing mass resolution (Just et al. 2011). The  \ndefinition 𝑏 min = 𝐺𝑀 BH / 𝑣 2 BH , combined with 𝑏 max = 𝜖 𝑔 , as in Tremmel et al. (2015), tends to sink the black holes too slowly.  \n(ii) The maximum unresolved impact parameter, 𝑏 max , unres , is independent of mass resolution and depends only on the gravitational softening and halo scale length, for which we provide a relation in Eq. (17). For most large-volume cosmological simulations and haloes hosting SMBHs, 𝑏 max , unres ≈ 0 . 5 𝜖 𝑔 . We emphasise that the exact value also depends on the softening spline employed in the simulations. Our results, in this case, are specific to the softening spline used in GADGET-2 (Springel 2005) and related codes.  \n(iii) We also obtain a relation between the number of particles required to accurately compute the local density and velocity dispersion as a function of the number of particles present in the halo (Eq. 19). We show that this substantially improves the error on these estimated quantities, compared to using a fixed number of neighbouring particles (e.g. 64 particles as in Tremmel et al. 2015). Directly counting the number of particles moving slower than the black hole (which are the ones contributing to dynamical friction) is more accurate than applying the Maxwellian approximation, as in Chen et al. (2022), which underestimates the density near the centres of cuspy Hernquist halos. As a result, the model of Chen et al. (2022) requires 𝑏 max , unres ≈ 6 𝜖 𝑔 to sink the black holes faster, although, at high resolution, this combination of parameters tends to sink black holes too fast.  \n(iv) We have quantified the effects of stochasticity on black hole orbits by varying the background halo mass resolution. We found that even in cases with 𝑀 BH / 𝑀 DM = 100 stochastic effects can alter a black hole's orbit, although the time it takes for the effects to become noticeable in terms of sinking time are longer. The magnitude of stochastic heating increases as the 𝑀 BH / 𝑀 DM ratio and 𝜖 𝑔 decrease. We have derived an order-of-magnitude correction for the effects of stochasticity, which, in effect, enhances dynamical friction (Eq. 20). This correction is, in theory, equivalent to increasing the dynamical massoftheblackholebythedarkmatterparticle mass. We emphasize that, due to the random nature of this problem, the recovered sinking timescales of black holes can become highly inaccurate, particularly for 𝑀 BH / 𝑀 DM ≤ 5. Our correction does, however, guarantee that the black hole sinks.  \n(v) We have validated our model on a number of test cases and levels of resolution, for circular and eccentric orbits. We showed that, for a set of random seeds, we on average accurately reproduce the black hole orbits in KETJU, with scatter in black hole sinking times increasing as 𝑀 BH / 𝑀 DM decreases.  \nOfcourse, our model has only been tested in spherical, equilibrium 𝜌 ∝ 𝑟 -1 halos and its validity may change when applied to different density distributions (particularly cored halos, see e.g. Antonini & Merritt 2012), flattened systems and systems undergoing a merger, which would introduce spatial and velocity anisotropy (Peñarrubia et al. 2004). On the other hand, a great advantage of this model and previous models of its kind, is that some dynamical friction is already captured in the simulations, and the correction we apply is only a local one. In other words, the accuracy of the recovered black hole motion should depend only on whether the Chandrasekhar formula represents well the impact of dynamical friction within ∼ 𝜖 𝑔 , which, depending on the spatial resolution of the simulation, may be a very small scale. Moreover, like models of Tremmel et al. (2015) and Chen et al. (2022), as well as repositioning schemes, our dynamical friction prescription does not conserve momentum. This may lead to unrealistic properties of the particles surrounding the black hole, particularly in regions where the dynamical friction correction is strong. We discuss the implications of this in detail in  \nAppendix C. Ultimately, testing the full range of validity of our model requires more sophisticated numerical setups and cosmological zoom simulations, as well as complementary high-resolution simulations with KETJU, all of which we will leave to future work.  \nCosmological hydrodynamics simulations are undoubtedly very powerful tools to study black holes and how they affect their host galaxies. Predictions can be made on black hole growth and merger rates, halo occupation distribution and populations of off-centre and wandering black holes. These predictions, however, have many modelling uncertainties. The treatments of black hole seeding, accretion, feedback, mergers, and interactions with other black holes all play an important role. In this work, we focused on the problem of modelling the effects of dynamical friction on black hole dynamics and introduced a subgrid model that has been calibrated on simulations where dynamical friction is fully resolved. We have shown that our model is able to accurately capture black hole sinking timescales when compared to KETJU, achieving higher performance compared to alternatives in the literature. It is therefore ready to be applied in the next generation of cosmological hydrodynamics simulations.", 'ACKNOWLEDGEMENTS': "The authors would like to thank Yueying Ni, Christian Partmann and Simon White for helpful discussions. This research was carried out using the High Performance Computing resources of the FREYA cluster at the Max Planck Computing and Data Facility (MPCDF, https://www.mpcdf.mpg.de) in Garching operated by the Max Planck Society (MPG). This project was developed as part of the Simons Collaboration on 'Learning the Universe'. This research benefited from the use of /e.pc/m.pc/c.pc/e.pc/e.pc (Foreman-Mackey et al. 2013), /n.pc/u.pc/m.pc/p.pc/y.pc (Harris et al. 2020), /s.pc/c.pc/i.pc/p.pc/y.pc (Virtanen et al. 2020) and /m.pc/a.pc/t.pc/p.pc/l.pc/o.pc/t.pc/l.pc/i.pc/b.pc (Hunter 2007).", 'DATA AVAILABILITY': 'The data analysed in this article can be made available upon reasonable request to the corresponding author.', 'REFERENCES': 'Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 319, 168', 'APPENDIX A: THE DIFFUSION COEFFICIENTS FOR A HERNQUIST HALO': 'In Fig. A1, we demonstrate how the dynamical friction forces compare to the stochastic forces for our fiducial Hernquist halo, with 𝑀 BH / 𝑀 DM = 1000. The diffusion coefficients were computed through explicitly evaluating Eqns. (10), (11) and (12), assuming no softening. It can be seen that the dynamical friction force due to the background particle mass (blue line) is of similar order of magnitude as the second-order diffusion coefficients (red and orange lines), which justifies our approximation of the stochastic force in Eqn. (20). Note that the second-order diffusion coefficient in the perpendicular direction, 𝐷 [( Δ 𝑣 ⊥ ) 2 ]/( 2 𝑣 ) , corresponds to 4 possible directions in which the stochastic force can act, while the parallel second-order diffusion coefficient, 𝐷 [( Δ 𝑣 ∥ ) 2 ]/( 2 𝑣 ) , acts only in one direction (the direction of motion of the black hole).  \nThe dashed red lines show how 𝐷 [( Δ 𝑣 ∥ ) 2 ]/( 2 𝑣 ) would change if the mass of the background particles were to increase (reducing  \n)  \n(  \nFigure A1. The magnitude of the diffusion coefficients for a 10 9 M ⊙ black hole on a circular orbit inside a 10 13 M ⊙ Hernquist halo, with the black holeto-background particle mass ratio 𝑀 BH / 𝑚 ∗ = 1000. The black line shows the magnitude of dynamical friction due to the mass of the black hole, the blue line shows additional dynamical friction arising from the mass of the background particles, 𝑚 ∗ . The red line gives the magnitude of the stochastic force in the direction parallel to the direction of motion of the black hole, and the orange line displays the magnitude of the stochastic force perpendicular to the direction of motion. The dashed red lines indicate how the magnitude of the second-order parallel diffusion coefficient would change if the background particle mass increases (and the mass ratio decreases).  \n<!-- image -->  \n𝑀 BH / 𝑀 DM ), increasing the impact of stochasticity relative to that of dynamical friction.', 'APPENDIX B: MAXIMUM UNRESOLVED IMPACT PARAMETER IN THE TREMMEL MODEL': 'As in Section 4.2, we here derive the optimal value of 𝜂 for the case 𝑏 min = 𝐺𝑀 BH / 𝑣 2 BH . This is shown in Fig. A2, resulting in a relation  \n𝜂 T15 = min [ 1 . 41 ( 𝜖 g / kpc ) -0 . 04 10 0 . 08 ln 2 ( 𝜖 g / kpc ) , 3 . 0 ] . (B1)  \nAgain, we emphasise that this relation is only applicable for the softening spline used in GADGET and AREPO.  \nInterestingly, the relation is clearly offset for the cases 𝑀 BH / 𝑀 DM = 50 and 100, unlike what we obtained under a different definition of 𝑏 min , suggesting different values of 𝜂 are required as the resolution changes. This is reflected in the left panel of Fig. 2, where the simulation with 𝑀 BH / 𝑀 DM = 50 requires stronger dynamical friction to reduce the sinking time. In the simulation with 𝑀 BH / 𝑀 DM = 100, the black hole sinking is also delayed early on, though it subsequently catches up at higher resolution levels due to the increased ellipticity of the orbit, induced by noise. As the local velocity dispersion appears to be an important component of 𝑏 min , we are able to achieve higher convergence across resolution levels with our method compared to Tremmel et al. (2015).', 'APPENDIX C: MOMENTUM CONSERVATION': 'Similarly to the models of Tremmel et al. (2015) and Chen et al. (2022), as well as repositioning schemes, our dynamical friction model does not conserve the momentum of the simulated system. This is demonstrated on the left panel of Fig. A3, where we show the fractional change of the magnitude of momentum compared to the initial conditions. It is evident that in the lower-resolution simulations (where the softening is large), the discrepancy is more significant and the change in momentum can reach ∼ 10 per cent of the initial value by the time the black hole sinks (at ∼ 5 Gyr). After the black hole reaches the centre of the halo, the change in the momentum of the system can become extremely high, especially at low resolution. This is because the black hole is not fully stable at the halo centre and, as long as 𝑏 min is smaller than 𝑏 max , unres = 𝜂𝜖 𝑔 , the deceleration of the black hole in high-density regions computed with our model is extremely large. For the high-resolution case, 𝑀 BH / 𝑀 DM = 5000 and 𝜖 𝑔 = 0 . 2 kpc, since the central regions are well resolved and 𝑏 min becomes greater than 𝑏 max , unres (setting the dynamical friction correction to 0) the black hole is able to stay at the centre of the halo and the momentum thereafter is essentially conserved.  \nAs the black hole (or any massive infalling object) is slowed down due to dynamical friction, it loses angular momentum, while the surrounding particles would, in turn, gain momentum. This process has previously been invoked as a mechanism for transforming cold dark matter cusps in the centres of haloes into constant-density cores (El-Zant et al. 2001, 2004; Cole et al. 2011). In the middle panel of Fig. A3, we demonstrate that our simulation with KETJU of a 10 9 M ⊙ black hole sinking to the centre of a 10 13 M ⊙ Hernquist halo (black solid and dashed lines) also predicts a shallower dark matter density profile once the black hole has sunk to the halo centre. An increase in the velocity dispersion is also evident below a radius of 800 pc (where the initial Hernquist profile contains 10 9 M ⊙ in dark matter). In coloured lines, we show how these halo properties differ when our dynamical friction model is applied. When the softening is small (< 400 pc), the flattening of the density profile is partially recovered (as momentum exchange still occurs on the scales where dynamical friction is resolved) and the increase in the velocity dispersion of the surrounding particles can also be seen as the surrounding particles are heated by the black hole. When the softening is large (> 400 pc), the slope of the density profile does not substantially change and neither does the velocity dispersion. Although the lower-resolution cases do have a central density core, this is primarily induced by the gravitational softening (Barnes 2012), while the background particle kinematics remain cold.  \nIn practice, the lack of momentum conservation would lead to underestimated velocity dispersion of background stars and dark matter in the vicinity of the black hole, with the effect being strongest when the gravitational softening is large, for very massive black holes and in dense regions, where the dynamical friction is strong. To circumvent this issue, one could, in principle, artificially inject momentum into the surrounding particles such that total momentum is manifestly conserved. It is, however, unclear how to distribute this momentum, such that the resultant dynamics of the surrounding particles faithfully reflect momentum transfer in dynamical friction. Perhaps, a weighting scheme for particles with distance 𝑑 < 𝜂𝜖 𝑔 of the black hole that follows Eqn. (9) of Ma et al. (2023) could be used, though this will likely break down at low resolution if the black hole is not stable at the centre of the halo (see the left panel of Fig. A3), leading to an unnecessarily large cumulative momentum injection into the surrounding matter. For now, it is clear that spatial and kinematic properties of particles surrounding a massive black  \nFigure A3. Left: fractional change in the magnitude of the total momentum of the simulated system (a 10 13 M ⊙ Hernquist halo with a 10 9 M ⊙ black hole) when our dynamical friction model is applied. Different colours correspond to different background particle resolution and softening (see Table 1). The black holes sink at ∼ 5 Gyr. Middle: the change in the density profile after 1 Gyr (solid line) and 8 Gyr of evolution (dashed line) for the 4 highest-resolution cases from the left panel. Black lines show the results of the simulation with KETJU (where momentum is essentially conserved, with Δ | 𝑝 | /| 𝑝 0 | ∼ 10 -8 ). The grey dotted line is the analytical Hernquist density profile. The grey shaded band covers the region of the halo containing 10 9 M ⊙ of dark matter in the initial conditions, equal to the black hole particle mass. The arrows show the gravitational softening used in each simulation. Right: as the middle panel, but showing the change in 3-dimensional velocity dispersion.  \n<!-- image -->  \nFigure A2. As Fig. 6, but using the Tremmel et al. (2015) definition of 𝑏 min .  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n]  \n1  \ns  \nm  \nk  \n[  \n)  \nr  \n(  \nhole in galaxy centres cannot be taken at face value in simulations where momentum exchange due to dynamical friction is not fully accounted for.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.'}