Dynamical friction of a massive black hole in a turbulent gaseous medium

The orbital decay of massive black holes in galaxies in the aftermath of mergers is at the heart of whether massive black holes successfully pair and merge, leading to emission of low-frequency gravitational waves. The role of dynamical friction sourced from the gas distribution has been uncertain because many analytical and numerical studies have either focussed on a homogeneous medium or have not reached resolutions below the scales relevant to the problem, namely the Bondi-Hoyle-Lyttleton radius. We performed numerical simulations of a massive black hole moving in a turbulent medium in order to study dynamical friction from turbulent gas. We find that the black hole slows down to the sound speed, rather than the turbulent speed, and that the orbital decay is well captured if the Bondi-Hoyle-Lyttleton radius is resolved with at least five resolution elements. We find that the larger the turbulent eddies, the larger the scatter in dynamical friction magnitude, because of the stochastic nature of the problem, and also because of the larger over- and under-densities encountered by the black hole along its trajectory. Compared to the classic solution in a homogeneous medium, the magnitude of the force depends more weakly on the Mach number, and dynamical friction is overall more efficient for high Mach numbers, but less efficient towards and at the transonic regime.


Introduction
The evolution of black holes (BHs), their ability to accrete surrounding matter, or to undergo mergers is directly linked to their dynamical evolution, their position in the galaxy, and their velocity relative to the ambient medium. In general, BHs move in the potential of their host galaxy. More locally, they exchange momentum and energy with the interstellar medium and the stars through dynamical friction.
Dynamical friction can be defined as the interaction of a massive object with its own gravitational wake. A massive object moving through lighter bodies or a gaseous medium attracts the surrounding matter, creating an over-density (wake) behind it. The wake in turn exerts a gravitational force on the massive object and slows it down. In this way, momentum and kinetic energy are transferred from the massive object to its over-dense wake.
Dynamical friction is an important force in many different contexts, including the evolution of BH binaries (Yu 2002;Dotti et al. 2007;Mayer 2013;Dosopoulou & Antonini 2017;Li et al. 2020b) and X-ray binaries (Iben & Livio 1993;Hurley et al. 2002), the orbital evolution of black holes in galaxies (Volonteri & Perna 2005;Bellovary et al. 2010;Tremmel et al. 2017), the orbital decay of a galaxy satellite in a galaxy cluster (Colpi et al. 1999;Fujii et al. 2006;Ogiya & Burkert 2016), among others. The detection of gravitational waves from supermassive BHs, which is possible with the upcoming Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017) and through pulsar timing arrays (Jenet et al. 2004(Jenet et al. , 2005, has made the orbital decay of supermassive BHs in a galaxy merger remnant a topic of special in-terest, as one needs to determine under which conditions supermassive BHs could coalesce and emit gravitational waves in less than a Hubble time. Furthermore, supermassive BHs can accrete efficiently and for a prolonged time from their surroundings, and exert feedback on the host galaxy, which is considered of paramount importance to regulate the baryon content in halos and galaxies, and to set up diverse galaxy properties (e.g. Granato et al. 2004;Di Matteo et al. 2005;Croton et al. 2006;Bower et al. 2006;Sijacki et al. 2007;Dubois et al. 2010Dubois et al. , 2012Dubois et al. , 2016Beckmann et al. 2017).
There is increasing theoretical support for wandering massive BHs within low-mass galaxies (Bellovary et al. 2019;Pfister et al. 2019;Lapiner et al. 2021;Ricarte et al. 2021), which might correspond to spatially offset active galactic nuclei as observed with radio emission in a number of dwarfs (Reines et al. 2020). This population of wandering BHs originate from galaxies tidally perturbed by a merger (Bellovary et al. 2021), can be scattered within galaxies by their multiphase structure (Pfister et al. 2019), or can be produced by gravitational recoil (Blecha et al. 2011). Therefore, accurately modelling the BH dynamics, which should consider the effect of dynamical friction, is key to BH mass growth and feedback (Lapiner et al. 2021;Bahé et al. 2022), their sinking into galaxies (Tremmel et al. 2017; Bartlett et al. 2021;Ma et al. 2021), and, hence, of their merger rate Li et al. 2020a; Barausse et al. 2020;Chen et al. 2022;Kunyang et al. 2022).
Dynamical friction was first described analytically by Chandrasekhar (1943) who developed the theory for a collisionless medium, such as stars or dark matter particles. Then Ostriker (1999) proposed an analytical solution for a uniform gaseous medium using time-dependent linear perturbation theory. Ostriker (1999) also showed, by comparing her work to the work of Chandrasekhar, that the gaseous dynamical friction is less efficient than the collisionless dynamical friction in the subsonic regime but much more efficient in the transonic regime where it sharply peaks. For large Mach number dynamical friction force is the same in a gaseous or collisionless medium. Recent insights have shown that dynamical friction from gas and stars both play a crucial role in allowing orbits of BH to decay in galaxies (Pfister et al. 2019;Chen et al. 2022), with gas dynamical friction particularly important at high redshift.
The non-linear evolution of dynamical friction has been studied using idealised numerical experiments. Early work was done by Ruffert & Arnett (1994) and Ruffert (1996) who studied the evolution of the wake and the resulting force on the BH using uniform background gas density and sound speed. Like most work in the field, these experiments were conducted in the frame of the BH, with the background gas moving at a fixed Mach number at all times. During this work, it became apparent that, as predicted in analytic models (Cowie 1977), wakes are unstable, a phenomenon further explored in Foglizzo et al. (2005). Since then, studies have investigated the impact of various physical phenomena on dynamical friction, such as for example density gradients (MacLeod & Ramirez-Ruiz 2015) or the impact of BH radiative feedback on wake structure (Park & Bogdanović 2017. As most of the force onto the BH is produced by density perturbations very close to the BH, resolving dynamical friction in large-scale cosmological simulations is numerically very expensive. Instead, it can be included in simulations using analytic sub-grid models (Dubois et al. 2013Tremmel et al. 2015;Pfister et al. 2019;Ni et al. 2022). While not without its challenges (Korol et al. 2016;Beckmann et al. 2018;Morton et al. 2021), this method has allowed for crucial improvements in reliably calculating the orbits of massive BHs in galaxies.
In this paper we study dynamical friction from a gaseous medium. The work presented here bridges the gap between isolated numerical experiments and full galaxy scale simulations to investigate the impact of small-scale turbulence on the dynamical friction force. The interstellar medium being turbulent, the conditions adopted here are more realistic than the thoroughly investigated, but idealised, case of a homogeneous gaseous medium.
The paper is organised as follows. In Sect. 2 we recall the theoretical background and provide definitions and modifications for our study. In Sect. 3 we introduce the numerical models. In Sect. 4 we describe our simulations and their initial conditions. Results are presented in Sect. 5. Section 6 is devoted to discussion and conclusions.

Theory of dynamical friction
The theory of dynamical friction in a gaseous medium was first based on linear theory under the assumption of a steady state. In the supersonic case, Dokuchaev (1964), Ruderman & Spiegel (1971), and Rephaeli & Salpeter (1980) obtained the following formula for the dynamical friction force exerted on a massive perturber: where G is the gravitational constant, M BH is the mass of the perturber, ρ 0 the background gas density, v rel the relative velocity of the perturber with respect to the background gas, and r min and r max are respectively the size of the perturber and the characteristic size of the surrounding medium. In many situations, r min and r max are poorly defined but due to the logarithm, large variations in the two characteristic sizes leads to a small difference in the overall force. In the subsonic case, Rephaeli & Salpeter (1980) argued that the dynamical friction force is null due to the symmetrical distribution of the gas around the massive object. However, this hypothesis shows a discontinuity between the supersonic and subsonic regimes, when the dynamical friction force, maximal just above Mach number suddenly goes to zero. This nonphysical result led Ostriker (1999) to review the problem under a time-dependent assumption.
The time-dependent theory allows to calculate a non-zero dynamical friction force in the subsonic regime: and to take into account the development of the wake in the supersonic regime: where M = v rel /c s is the Mach number, and c s is the sound speed. The analytic force is therefore discontinuous at the sonic point, but a linear interpolation can be used to link the two regimes. According to Ostriker (1999), the characteristic size of the medium r max is equal to the size of the wake, r wake . This, in turn can be calculated from v rel t (where t is the time after which the perturber has started to exert its force on the background medium) under the assumptions that (v rel − c s )t is bigger than the perturber size r min and that (v rel + c s )t is smaller than the surrounding medium size r max .
For the turbulent medium we are interested in, some additional definitions are needed. In this paper, we use two different definitions of the Mach number. Firstly the classic Mach number M = v rel /c s . Throughout the text, categorisations such as 'supersonic' (M > 1), 'transonic' (M ∼ 1), and 'subsonic' (M < 1) refer to the classic Mach number. Secondly a turbulent Mach number M turb = v rel /v eff where v eff = c 2 s + v 2 rms is an effective gas velocity that captures the contribution of the gas sound speed and of the turbulent velocity (v rms is the root mean square velocity of the gas).
For a turbulent medium, we expect the time it takes for the stirring of the turbulence to erase the wake over-density to determine the typical length of the wake, r wake . The turbulence will erase the wake on a characteristic timescale t cross = r BHL /v eff , which is equal to the time taken by the effective turbulent velocity v eff to cross the width of the wake, taken equal to the Bondi-Hoyle-Lyttleton (BHL) radius r BHL , defined as: Hence, the length of the wake is expected to simply be r wake = t cross v rel = r BHL M turb when fully established. Taking early times into account, we have:

Magneto-hydrodynamics simulations with forcing turbulence
The simulations presented here were performed with the ramses code (Teyssier 2002), solving for ideal magnetohydrodynamics (MHD; Fromang et al. 2006) with self-gravity. The gas evolution was obtained with a MUSCL-Hancock scheme which used a second-order Godunov method with a minmod limiter on the linear reconstruction of the conserved quantities at cell interface. The MHD flux at each cell interface was obtained with the approximate HLLD (Harten-Lax-Van Leer-Discontinuities) Riemann solver (Miyoshi & Kusano 2005). The induction equation that evolves the magnetic B field was solved with the constrained transport algorithm that guarantees ∇.B = 0 at machine accuracy (Evans & Hawley 1988;Teyssier et al. 2006). The massive perturber was evolved with a particle-mesh method using a cloud-in-cell interpolation. The total gravitational potential, also accounting for the gas contribution, was obtained with a multi-grid solver for the Poisson equation (Guillet & Teyssier 2011). Simulations presented here were run using a uniform grid, as the turbulent nature of the gas makes numerical savings from using adaptive mesh refinement very small. In addition, a uniform grid avoids errors in the Poisson solver, that can generate a non-negligible self-force at coarse-to-fine boundaries near the BH (Bleuler & Teyssier 2014;Zhu & Gnedin 2021) due to its strong gravitational potential.
We used a specific model to drive turbulence in the interstellar medium, based on the implementation done by Commerçon et al. (2019) with further details of the implementation documented in Schmidt et al. (2009) andFederrath et al. (2010). Kinetic energy was injected using the Ornstein-Uhlenbeck process (Eswaran & Pope 1988), which creates a turbulent forcing generated in Fourier space, using a combination of compressive (1/3) and solenoidal (2/3) modes. This forcing was defined by the wavenumber k turb , the power spectrum of shape 1 − (k − k turb ) 2 , and an ad hoc boost factor of the driving force, which sets the amplitude of the turbulent velocity. In the following set of simulations we set the boost factor to a constant value for any choice of central k turb , which turns into a smaller rms velocity for larger k turb . The typical range of values of rms velocities experienced in these simulation were within v rmx = 2−4.6 km s −1 (for k turb = 16 to 4 respectively). The turbulence wavenumber k turb , referred to as k in the following, is a fraction of the box size and represents the mean size of turbulent eddies. The auto-correlation time was set to 0.6 Myr.
The gas could cool down to temperatures as low as 10 K using the cooling function L(ρ, T ) from Sutherland & Dopita (1993) for temperatures above 10 4 K and the cooling due to carbon, oxygen, and hydrogen for temperatures below 10 4 K. We assumed solar metallicity and a typical Milky Way UV flux for ionisation counts. The gas was a perfect gas with an adiabatic index γ = 5/3 and a mean molecular weight µ = 1.4. We neglected any gas accretion onto the BH, and hence its feedback, as well as star formation and feedback from stars occurring in the interstellar medium.

General setup
In order to study the effect of the multi-phase gas structure seeded by the turbulence we have chosen to set an initial gas density of ρ 0 = 10 cm −3 , within the typical range of the neutral medium so that the gas got thermally unstable and formed structures. The turbulence energy injection provided an additional heating leading to a temperature equilibrium about T 350 K with a standard deviation of σ T 345 K and a sound speed about c s 1.2 km s −1 with a standard deviation of σ c s 0.6 km s −1 . The magnetic field at equilibrium was about B 2.4 µG, corresponding to an Alfvén velocity about v A 1.4 km s −1 . In our simulations, the Alfvénic Mach number (M A = v rms /v A ) remained above one, that is in the supersonic or nearly transonic regime. Although magnetic fields can play a significant role in the anisotropic shape of the turbulence (e.g. Goldreich & Sridhar 1995;Cho & Vishniac 2000;Beattie & Federrath 2020), we deferred the investigation of its effect on dynamical friction to future work.
In idealised simulations, the BH mass and velocity cannot be chosen arbitrarily as there are a series of constraints that must be fulfilled to combine sufficient resolution with a reasonable computational time, within a physically interesting set-up.
The resolution was set by the requirement to resolve the BHL radius r BHL . Resolving r BHL was necessary to capture the gas overdensity formed by the BH gravitational potential, and, hence, the dynamical friction force exerted by the gas onto the BH (Beckmann et al. 2018, and references therein).
Furthermore, we had to take into account that r BHL increases over time as the BH slows down. A maximum value is reached for r BHL when the BH velocity is similar to the effective gas velocity. To study the problem presented here, r BHL must be kept at all times between a few times larger than the resolution and a few times smaller than the box size. Figure 1 shows the value of r BHL for our reference case for the gas properties and for different values of the BH mass and BH relative velocity. Based on this, the initial BH velocity was set to about 3 times that of the gas. This allowed us to study a range of Mach numbers while avoiding large variations of r BHL over time.
Given the typical v eff of the gas distribution, this series of constraints led to 10 4 M for the BH mass and 40 pc for the box size as acceptable values. To avoid excessive computational times we had also calculated the expected theoretical dynamical friction time, and we had checked the BH-to-gas mass ratio and the Jeans mass to avoid gas collapse, but these were less restrictive parameters. For the Jeans mass calculation we noted that the relevant velocity to consider is the effective velocity v eff ≥ 3.3 km s −1 , which gives a Jeans mass of about 10 5 M , larger than the total gas mass in the box (2 × 10 4 M ). The injected turbulence was therefore efficient enough to prevent gas collapse; this was confirmed by the simulations where we have not observed any (complete) gas collapse. Small transient gas over-densities can appear, but they are rapidly dispersed by the turbulence. About 10 Myr (corresponding to a turbulence box crossing time) are needed for the turbulence to settle and for the mean gas parameters, temperature, pressure, to stabilize. After this first period of 10 Myr, a BH was added, with a mass of 10 4 M and an initial velocity of 15 km s −1 . Periodic boundary conditions allowed the BH (and the gas) to cross the box many times throughout one of our simulations. To avoid the BH crossing its own wake, the BH had an initial trajectory with an angle of 30 • with respect to the x-axis and 60 • with respect to the z-axis.
In this paper, we study the impact of resolution, turbulent wavenumber, and stochastic variation on the BH trajectory over time. The fixed parameters used for all simulations are summarised in Table 1, while parameter choices specific to a given simulation are reported in Table 2.

Run0
Run0 will serve as a reference case throughout this paper. It has been selected to have an intermediate level of turbulence, with v rms = 3.2 km s −1 , a wavenumber of k = 8, and r BHL that is resolved throughout the simulation. As discussed in Sect. 4.1, the BHL radius increases with time as the BH velocity decreases. With a resolution of ∆x = 0.16 pc for Run0, r BHL evolves from 2.4∆x in the supersonic regime at the beginning of the simulation to 26∆x in the subsonic regime at the end of the simulation.
A projection of the gas density at time t = 10 Myr can be seen in the central panel of Fig. 2. For this simulation, typical eddy sizes are about 5 pc, which is slightly larger than r BHL , which ranges from 0.4 to 4 pc.

Turbulent seed
Turbulence was initialised with a random seed, TSeed, which determines the exact instance of density peaks and troughs in the initial conditions of the simulation. For simulations with different TSeed, the gas has broadly the same appearance, with over-and under-densities having on average the same size and amplitude, but the local conditions encountered by the BH as it crosses the box differ. To test the statistical variation of results, we conducted a set of four simulations that had an identical setup to Run0 but were initialised with a different random seed. These simulations are called TSeedA, TSeedB, TSeedC, and TSeedD. The aim of this series of simulations is to estimate the effect of the stochastic processes at work and to assess the difference between the general trends of the gas and BH evolution and any random event that could happen.

Resolution
To explore the response of the system and the dynamical friction force with respect to how well the BHL radius is resolved, we had set up simulations for a range of resolutions, varying between 0.08 ≤ ∆x/pc ≤ 0.62: -with a ∆x = 0.62 pc resolution (Res0.62pc), the BHL radius is unresolved (r BHL ≤ ∆x), -with a ∆x = 0.31 pc resolution (Res0.31pc), the BHL radius is marginally resolved (r BHL ∆x) at the beginning of the run and resolved (r BHL ≥ 2.5∆x) after 120 Myr, -with a ∆x = 0.16 pc resolution (Run0), the BHL radius is resolved, -with a ∆x = 0.08 pc resolution (Res0.08pc), the BHL radius is resolved, with r BHL equal at least to 5∆x.

Turbulence wavenumber
To assess the impact of different sizes of turbulent eddies, we varied the turbulence wavenumber, that is the size of the individual turbulent regions as a function of the box size. Values tested here ranged from k = 4 to k = 16, leading to a driving scale of the turbulence from 10 pc to 2.5 pc: -with a k = 4 turbulent wavenumber (k4, k4_TSeedA, k4_TSeedB simulations), the driving scale of the turbulence is about 10 pc, that is at least 3 times the BHL radius, -with a k = 8 turbulent wavenumber (Run0, TSeedA, TSeedB, TSeedC, TSeedD, but also the simulations with various resolution Res0.08pc, Res0.31pc, Res0.62pc), the driving scale of the turbulence is about 5 pc, that is always larger than the BHL radius but ranging from 13 times larger at the beginning of the simulations to the same order of magnitude at the end, -with a k = 16 turbulent wavenumber (k16, k16_TSeedA, k16_TSeedB simulations), the driving scale of the turbulence is about 2.5 pc, that is the order of magnitude of the BHL radius. Snapshots of the resulting density distribution of simulations with k = 4 and 16 can be seen in Fig. 2 (respectively the left and right panels).

BH deceleration
When placed in a turbulent background medium, an over-dense wake develops downstream of a moving BH (see Fig. 3), which creates dynamical friction that is able to slow an originally supersonic BH (with initial Mach number M ini = 12.4; this initial Mach number is derived from the average sound speed calculated from the simulation at the first output after BH insertion, at t = 15 Myr) down to the transonic regime (see Fig. 4, left panel) where the evolution stalls. This is consistent with expectations from analytical theory, where dynamical friction in the subsonic regime becomes negligibly efficient. We can note that the BH slows down to the sound speed, rather than to the turbulent (rms) velocity of the gas (here and after the rms velocity is measured for the gas outside a sphere of 5 pc radius centred on the BH to avoid including local increase in velocity caused by A217, page 4 of 14 the growing overdensity of gas around the BH at late times in transonic regime). When the BH motion becomes transonic with respect to the turbulent Mach number (i.e. when M turb 1), there is no noticeable effect.
During the simulation, the gaseous wake grows from a long and extended structure during the early supersonic phase to an almost spherical overdensity around M ≈ 1. Figure 3 shows the wake evolution after t = 30, 85, 120 Myr (supersonic regime), and 265 Myr (transonic regime) respectively.
The BH motion follows three main phases: -Phase 1. During the first 110 Myr (from the BH injection at 10 Myr to about 120 Myr), the BH is already slowing down but at a low rate of about 0.04 km s −1 Myr −1 (Fig. 4, right panel, phase highlighted in red). During this phase, the BHgas interaction is at work but appears to be limited and the dynamical friction not fully captured. This trend is confirmed with the resolution study in Sect. 5.3. -Phase 2. During the next 100 Myr (from 120 Myr to 220 Myr), the decrease in BH velocity more than doubles in comparison to the first phase with a rate of 0.088 km s −1 Myr −1 (Fig. 4, right panel, yellow). The BH-gas interaction is reinforced and the dynamical friction is now well captured.
-Phase 3. After 220 Myr, the BH enters the transonic regime. During this phase, the BH velocity oscillates around the value of the sound speed but it does not decrease any more. At that time dynamical friction becomes inefficient and no further trends are observed until the end of the simulation. The BH velocity oscillates around 1.6 km s −1 (Fig. 4, right panel, green), close to the gas sound speed of about 1.14 km s −1 (while the rms velocity is larger at about 3 km s −1 ). We discuss the trends in detail in the parameter study sections, as trends can be quantified and explained more robustly when comparing Run0 to other simulations. Here we briefly note that the first two phases are related to how well the Bondi-Hoyle-Lyttleton radius is resolved. Early on, when r BHL is barely resolved, with less than ∼5 resolution elements, dynamical friction is not fully captured (first phase). Later on, r BHL /∆x increases and dynamical friction becomes more efficient.

Development of the wake
In this section, we focus on the wake developing behind the BH, which is the over-dense gas responsible for the BH deceleration. In order to quantify the wake evolution, we measure the average gas density in a cylinder of radius 2 pc and length 20 pc just behind the BH, aligned with the BH velocity vector. The BH is injected at 10 Myr and the first available results shown here are at 15 Myr. The free fall time of gas in the BH potential for a BH of 10 4 M travelling at 15 km s −1 in a gas with an effective velocity of 3 km s −1 is about 0.04 Myr. By employing a sampling frequency of 5 Myr, we are therefore not following the set-up of the wake but its further development throughout the simulation. To have a reference value before the BH injection, we also measure the density in the same cylinder positioned according to the BH initial position and velocity direction in the first output of the simulation. Figure 5 shows the density evolution of the wake compared to the BH and gas velocity evolution for the Run0 simulation. The black dashed vertical lines are for the first BH output, 5 Myr after its injection (t ∼ 15 Myr) and the last output in the supersonic regime respectively (t ∼ 220 Myr). Before the BH injection, the density measured in the cylinder is about 10 cm −3 , in agreement with the gas mean density in the full box. After the BH injection the density quickly doubles to about 20 cm −3 and continues to increase from there. The wake evolution is very different in the supersonic and in the transonic regimes: in the supersonic regime the wake density grows slowly at a rate of 2.6 cm −3 per 10 Myr, whereas, in the transonic regime, the wake density shows large oscillations correlated with the BH relative velocity. Peaks in density correspond to low BH velocity and vice versa. Despite the oscillations, the average density of the wake continues to increase due to the gravity of the BH.

Sensitivity to initial turbulence seed
A turbulent gaseous medium is by definition a chaotic system, particularly sensitive to initial conditions. In order to assess the reliability and the stability of our results we run a set of simulations with the same physical and numerical conditions as Run0 but changing the turbulence seed. These are the simula-tions named TSeedA, TSeedB, TSeedC, and TSeedD for which the evolution of the BH relative velocity is presented in Fig. 6.
The results are quite homogeneous. TSeedA and TSeedB show a very similar evolution with a transition from phase 1 to phase 2 a little less pronounced than in Run0 but a transition to phase 3 at the same time (220 Myr). TSeedC and TSeedD on the other hand have a more noticeable transition between phase 1 and phase 2 leading to a slightly earlier transition to phase 3, about 10 to 20 Myr earlier (see vertical lines in Fig. 6). These differences are quite small and the global evolution of the BH due to dynamical friction is robust against small stochastic changes.
The consistency of our results was confirmed by repeating this analysis for simulations with other turbulence wavenumbers k = 4 and k = 16, which also showed good consistency. These results are available in Appendix A.

Importance of resolution
The resolution of our main simulation Run0 (0.16 pc) was chosen in order to resolve the BHL radius and indeed we were able to capture dynamical friction. To assess the effect of resolution, we here compare Run0 to one simulation with a higher resolution, one simulation with a lower resolution and one where the BHL radius is not resolved at any time. These are the simulations named Res0.08pc, Res0.31pc, and Res0.62pc respectively. The evolution of the BH relative velocity and the ratio between the BHL radius and resolution for this set of simulations is presented in the left and right panels of Fig. 7 respectively.
As can be seen in Fig. 7, the resolution has a significant impact on the results. When ∆x is comparable to or larger than the BHL radius, dynamical friction is not captured and the BH experiences little deceleration. This is the case for the Res0.62pc simulation, which has r BHL ≤ 1.2∆x at all times, and is unable to build an overdensity within r BHL . And this is also the case for the first part of the Res0.31pc simulation, until 130 Myr, with a BHL radius lower than 1.5∆x. When r BHL < 1.5∆x (Fig. 7, right A217, page 6 of 14 Fig. 4. Run0 velocities evolution. Left: evolution of the BH relative velocity (blue), the gas effective velocity (red, excluding a sphere of radius 5 pc around the BH), and sound speed (yellow) over time for Run0. Right: evolution of the BH relative velocity with superimposed fits to its time evolution. The dashed vertical lines in both panels mark the three phases of the BH deceleration. In the first phase the BHL radius is not well resolved, in the second the BHL radius is well resolved and deceleration is stronger, in the third the BH is in the transonic regime and does not decelerate further. Fig. 5. Wake density evolution compared to the BH and gas velocity evolution for the Run0 simulation. The gas density is measured in a cylinder of 2 pc radius and 20 pc length just behind the BH, oriented along the BH velocity vector. In the measure of the gas effective velocity the gas within a sphere of radius 5 pc around the BH is excluded. The black dashed lines show the first BH output, 5 Myr after its injection and the last output in the supersonic regime respectively. In the supersonic regime density in the wake grows slowly and steadily, whereas, in the transonic regime, the density keeps increasing, but with large oscillations. panel horizontal black dashed line) dynamical friction is not effective and the BH deceleration is minimal. The turning point at 130 Myr for the Res0.31pc simulation, when dynamical friction becomes efficient, can clearly be seen in Fig. 7 (left panel, yellow curve).
When the BHL radius is larger than 1.5∆x, dynamical friction is active but significantly lower than when r BHL > 4.5∆x ( Fig. 7 right panel, horizontal black dotted line), at which point a (second) increase in slope can be seen in the evolution of velocities (vertical coloured dotted lines, left hand panel). We have already commented on this in Run0 with a change in slope of the BH velocity at 120 Myr (vertical red dotted line in Fig. 7), at which point the deceleration doubles to 0.088 km s −1 Myr −1 . The same happens in the Res0.31pc simulation at 230 Myr (vertical yellow dotted line in Fig. 7) with the same deceleration rate around 0.088 km s −1 Myr −1 after this time. This result is again confirmed in Res0.08pc, which always has r BHL > 5∆x, and where the deceleration rate is consistently ∼0.089 km s −1 Myr −1 , except during the first 50 Myr when the wake is first forming. We therefore conclude that a minimum resolution of r BHL /∆x > 5 is Fig. 6. Evolution with time of the BH relative velocity for different turbulence seed with a turbulence wavenumber k = 8. The vertical dashed lines mark the transition from super-to transonic; simulations Run0, TSeedA and TSeedB transition exactly at the same time, therefore only one line is visible (red dashed). This shows that our results are robust to small stochastic changes in the initial conditions. required to correctly capture the magnitude of dynamical friction in a turbulent medium, with results converged for the range of resolutions probed here out to r BHL /∆x ≈ 50.
By studying the evolution of wake density with simulation time for simulations with different resolutions, we can show that the wake grows faster with higher resolution in the supersonic regime (see Fig. 8). However, when comparing the wake evolution as a function of Mach number, rather than simulation time, a global pattern emerges. To compare the density across simulations in a time-independent manner, we highlight in Fig. 8 the time at which the BH first becomes trans-sonic using vertical dotted lines. Comparing the wake density at these points in time shows that density in the wake is similar for the three different resolution runs when the BH reaches M = 1. In summary, in the supersonic regime, when increasing the resolution while keeping the same turbulence wavenumber, the density in the wake reaches the same level in a shorter time. In the transonic regime the mean peak density, ignoring the highest peaks and deepest troughs, appears to increase with resolution although in the Res0.31pc case we can not conclude whether the wake is still growing at the end of the simulation. The trend with resolution is consistent with gas being more easily bound to the BH A217, page 7 of 14 The dashed and dotted horizontal lines in the right panel represent respectively BHL radius equal to r BHL = 1.5∆x and r BHL = 4.5∆x. In the simulation where the BHL radius is unresolved (Res0.62pc), BH deceleration is very inefficient. In Res0.31pc, where the BHL radius is resolved after 120 Myr, we see at the same time that the BH starts to decelerate faster, with a second increase in deceleration occurring when r BHL = 4.5 at t = 240 Myr. As the resolution increases and the BHL radius is resolved with more elements (Res0.16pc and Res0.08pc) the effect of dynamical friction becomes more pronounced. at high resolution. In the transonic regime, the wake becomes a more spherical overdensity, which continues to grow in mass, although it never reaches the state of a stable sphere.

Influence of turbulence wavenumber
Depending on the surrounding medium, BH may experience turbulence on different scales. We focus here on the influence of this turbulence scale on the dynamical friction efficiency. For high turbulence wavenumber k, the BH passes through many density peaks but each passage lasts a short time. With low k the BH crosses a smaller number of bigger over-densities. The impact of changing wavenumber can be seen visually in Fig. 2. Figure 9 shows the BH and gas rms velocity evolution for three different turbulence wavenumbers, with each wavenumber represented by three or four different turbulent seeds. As the simulations are done with a constant turbulence rms forcing, they have somewhat different initial rms velocities depending on the wavenumber: about 2 km s −1 with k = 16, 3 km s −1 with k = 8, and 4.5 km s −1 with k = 4. In the calculation of the rms velocity a sphere of 5 pc radius has been excluded, in order to remove the contribution from gas infalling towards the black hole, which becomes large when the black hole reaches the transonic regime.
In the cases k = 4 and k = 8, the rms velocity remains roughly constant; whereas, in the case k = 16, it shows spikes up to 3.5 km s −1 . This is because the overdensity around the black hole becomes larger than 5 pc in this case, since the lesser effect of turbulence allows for the overdensity around the black hole to grow more easily.
The BH deceleration is slower in the k = 4 case than in the k = 8 and k = 16 cases because the wake is less stirred by the turbulence for high k, and can therefore persist for longer while growing somewhat more strongly. We speculate that this reduction in force for k = 4 occurs because a higher initial rms velocity makes the transfer of momentum and energy from the BH to the gas less efficient.
Performing the wake analysis on the simulations with different turbulence wavenumber, we show in Fig. 10 that the overdensity in the wake grows more strongly when the turbulence wavenumber is higher. The density reached by the wake at the end of the supersonic phase depends mainly on the turbulence wavenumber, and increases with it: from about n = 35 cm −3 for k = 4, via 70 cm −3 for k = 8 (Run0, Res0.08pc, and Res0.31pc) to about 150 cm −3 for k = 16. This confirms that larger eddies are more efficient in dissipating the over-density behind the BH. In summary, in the supersonic regime, when increasing the turbulence wavenumber while keeping the same resolution, the density in the wake reaches a higher level, that is, it is stronger, over the same period of time.
Once the BH reaches the transonic regime, the wake density evolution is rather chaotic with strong oscillations, as can be seen in Fig. 10. In general, although the BH does not slow down any more, the wake density keeps increasing (except in the k = 4 case). This increase is actually stronger than in the supersonic regime (equal in the k = 4 cases) and is weaker with a lower turbulence wavenumber, since the growing overdensity around the BH is denser, larger, and more stable with higher k.

Drag force as a function of Mach number
The last term in Eq. (3) includes ln(r max /r min ), which is sometimes referred to as the Coulomb logarithm. This term contains A217, page 8 of 14 Fig. 9. Evolution with time of the BH relative velocity (left panel) and the gas rms velocity (right panel, excluding a sphere of radius 5 pc around the BH) for 3 different turbulent wavenumbers. Different linestyles of the same colour probe the stochastic variation of our results. Red curves show results for simulations k4, k4_TSeedA, and k4_TSeedB, blue curves show Run0, TSeedA, TSeedB, TSeedC, and TSeedD and yellow curves show k16, k16_TSeedA, and k16_TSeedB. Deceleration is slower at low wavenumbers, where, in our simulations, the rms velocity is higher, and plausibly transfer of kinetic energy from the BH to the gas less effective. Results are robust to stochastic variation of the turbulent density field. Fig. 10. Wake density as a function of time for turbulent wavenumbers k = 4, k = 8 (Run0), and k = 16. The dotted vertical lines mark the transition from supersonic to subsonic regime in each case. The density in the wake grows faster and to higher values at high wavenumbers, where the gas is closer to homogeneity. At low wavenumbers the larger eddies in the turbulence and the higher rms turbulence erase the wake more easily. the only free parameters of the model. In this section, we compare the dynamical friction force experienced by the BH in our simulation to Eq. (3) for a range of different values of ln(r max /r min ). For easier comparison, we define the dimensionless drag force f (M) as where F DF is the dynamical friction force experienced by the BH. For the analytic values, F DF = F sup from Eq. (3), so To calculate F DF from our simulations, we compute the instantaneous acceleration of the BH by dividing its net velocity change during one simulation timestep by the length of that timestep.
The distribution of f (M) as a function of Mach number can be seen in Fig. 11. To compare this to the analytic function, we plot curves computed using Eq. (3) for different values of ln(r max /r min ). Fitting the values of ln(r max /r min ) is not trivial because of the intrinsic scatter in the value of f (M) introduced by the turbulence, so we prefer to bracket the distribution of data points from our simulation. We also add curves fitted on results for M > 4, keeping in mind that these fits are quite poor.
Scatter in the measurement of f (M) is introduced in two different ways. Firstly, the efficiency of the drag force depends on the background density, and therefore varies even at fixed Mach number as the BH traverses underdense or overdense regions. Secondly, the relative velocity between the BH and the background gas varies because the BH enters regions of higher or lower gas relative velocity, which in turn changes the Mach number we measure for our BH. This effects becomes stronger at lower Mach number, where velocity oscillations are more pronounced.
As can be seen in Fig. 11, the shape of f (M) is significantly flattened in comparison to the analytic value at constant ln(r max /r min ). As a result, the value of ln(r max /r min ) that allows the analytic values from Ostriker (1999;Eq. (3)) to best match our simulations decreases significantly with Mach number. Phrased differently, if we were to fit f (M) in a turbulent medium with Eq. (3) and a single ln(r max /r min ), the fit would under-predict f (M) at high Mach number and over-predict it at low Mach number. We explore this phenomenon further in Sect. 5.5.2.
The turbulent wavenumber k also affects the dispersion of the results: we observe a smoother evolution of f (M) with high k (see Fig. 11, bottom right panel for the k = 16 case) while the scatter increases with low k (see Fig. 11, bottom left panel for the k = 4 case). This can be quantified using the values of ln(r max /r min ) which envelope all points, as plotted as dotted (minimum ln(r max /r min )) and dashed (maximum ln(r max /r min )) lines in Fig. 11. For k = 4, ln(r max /r min ) ranges from 0 to 15, while the range tightens to 0.05 to 12 in the k = 8 case, and even further to 0.1 to 8 in the k = 16 case. Also shown on the figure are fitted curves for all values of f (M) for M > 4 (dashed-dotted lines). These fits do not reveal any strong dependence on k as the fit is generally poor due to the stochastic nature of the problem and the decrease in efficiency of dynamical friction with decreasing Mach number in our simulation.
When looking at the impact of resolution (Fig. 11, top left panel), dynamical friction becomes more efficient at higher resolution, as can be seen by the fact that the values of f (M) for ∆x = 0.31 pc (red) generally lie below those for ∆x = 0.08 pc (green). However, higher resolution simulations also show more scatter, and therefore require a broader range of ln(r max /r min ) to A217, page 9 of 14 envelop the data (0.05 < ln(r max /r min ) < 17 if ∆x = 0.08 pc, compared to 0.05 < ln(r max /r min ) < 10 if ∆x = 0.31 pc).
In conclusion, the dynamical friction forces as a function of Mach number on the BH in a turbulent medium have an intrinsically flatter shape than in the homogeneous case derived in Ostriker (1999). In a turbulent medium, at fixed M, the magnitude of the force increases for higher resolution and higher k.

Efficiency of the drag force in a turbulent medium
Given that our simulation has a uniform resolution, r min = ∆x is well defined for our simulations. This leaves r max as the only free parameter in ln(r max /r min ) to quantify the magnitude of the dynamical friction. We speculated as to the expected values of r max for a turbulent medium in Sect. 2. In this section, we fit the magnitude of the drag force measured from the simulation using Eq. (3) to measure r max , assuming r min = ∆x. We note that r max discussed here is not the physical extent of the wake behind the BH, which we shall refer to as r wake instead, but reflects the net deceleration of the BH from both the wake and any other gravitational forces on the BH due to turbulent over-densities.
The results for this fit as a function of Mach number can be seen in Fig. 12. Two things are immediately apparent: r max is strongly a function of Mach number, and the values found here range from unphysically small (r max ≤ r min ) to unreasonably large (r max l box ). The physical extent of the wake, r wake is difficult to measure in simulations due to the density structures inherent in turbulent gas, which make it difficult to differentiate between the overdensity caused by the BH and an overdensity that is simply part of the turbulent gas. However, it certainly cannot exceed the size of the box. This sets a physical upper limit on r wake of around 40 pc. The r max required to explain the force according to Ostriker (1999) is also clearly higher than the theoretical estimate for r wake in Eq. (5), which remains 2−3 pc throughout the simulation.
The magnitude of r max is a strongly increasing function of Mach number, which reflects the flattened distribution of f (M) as a function of Mach number in comparison to the homogeneous case (see Sect. 5.5.1 for a discussion). This is in clear contradiction with both the original prediction by Ostriker (1999), who predict that r max ∝ Mt, and the trend predicted in Sect. (2) where r max ∝ M −1 (see Eq. (5)). Extrapolating from our turbulent simulations to a smooth case (i.e. one for which k → ∞) is not trivial, as discussed in Appendix B. Nevertheless, using a simulation as close to Run0 as feasible, we confirmed that the k = 16 case presented here shows good convergence with a fully smooth case for M > 10, while at lower Mach numbers the turbulent case is always below the smooth run, highlighting that the suppression of the force at low Mach numbers is caused by turbulence (see Appendix B for details). Another value for comparison for a locally smooth medium can be taken from Chapon et al. (2011), who used the dynamical evolution of a transonic BH in a smooth galaxy in the range M = 1.3 to 0.8 to predict r max /r min ∼ 25, that is r max ≈ 4 pc for our resolution. This is in good agreement with our results for both high turbulent Mach numbers (see Fig. 12) and the smooth simulations (see Appendix B).
We conclude that the net long-term dynamical friction in a turbulent medium, once fully resolved, is somewhat more efficient than in the analytical case for strongly supersonic BHs (of M > 5) but significantly less efficient as the BH approaches the transonic regime. Finally, we note that the difference in r max /r min is exacerbated by the exponential in the functional form, and also that in the scale-free formalism of Ostriker (1999) r min is not fully defined, while in simulations r min is necessarily defined to be the resolution. Compared instead to a homogeneous case simulation as described in the Appendix B, dynamical friction in a turbulent medium, once fully resolved, is similar at high Mach numbers and less efficient close to and at the transonic regime.
As can be seen in Fig. 12 (top right and bottom panels), lower turbulent wavenumbers produce significantly more scatter in r max but no significant difference in the average r max at a given Mach number. This is in good agreement with the fact that the BH velocity evolution in Fig. 9 shows a small delay in deceleration for k = 4 and little difference between k = 8 and k = 16.
One caveat to this conclusion is that the magnitude of r max calculated using Eq. (3) is very sensitive to the values of the constants, such as the BH mass. ρ 0 is particularly difficult to define. Here it is taken to be the average background density, but obviously the BH passes a series of over and under-densities so ρ 0 is poorly defined in this context. Another hint that a turbulent medium cannot be well approximated by ρ 0 comes from the fact that there is a clear dependence of the maximum r max on the turbulent wavenumber k: higher wavenumbers lead to lower max-imum values of r max and less scatter (Fig. 12, bottom panels). This suggests that v rms also plays a role in determining the magnitude of the force. Overall, we conclude that dynamical friction due to a turbulent medium for strongly supersonic BH is more efficient than in the homogeneous case, but less efficient as the BH approaches the transonic regime.

Conclusions
In this paper, we presented a study of dynamical friction in a turbulent gaseous medium. We investigated how the deceleration of an initially supersonic BH depends on the properties of the local turbulent medium, thereby bridging the gap between isolated numerical experiments and full-scale galaxy evolution simulations. We find that: -In a turbulent medium, an over-dense wake develops downstream of the BH, which slows an originally supersonic BH (here, M ini ∼ 10) down to the transonic regime, where M = 1. -The efficiency of dynamical friction depends on how well the Bondi-Hoyle-Lyttleton radius is resolved in the simulations. For 1 < r BHL /∆x < 5, the dynamical friction is present but reduced in comparison to higher values. At higher resolution, the magnitude of the force is converged for the full range of resolutions tested here, 5 < r BHL /∆x < 50. -Even in a turbulent medium, evolution of dynamical friction is determined by the classic thermal Mach number M = v rel /c s rather than the turbulent Mach number M turb = v rel / c 2 s + v 2 rms which takes the turbulent rms velocity into account. While nothing special happens when M turb = 1, dynamical friction becomes very inefficient when M = 1 and the BH velocity evolution stalls at this point with no further slow-down observed in our simulations. -Dynamical friction is more efficient for higher wavenumbers of the background turbulence. In our simulations, the BH reaches the transonic regime after approximately 200 Myr for k = 16, rather than the 300 Myr required when k = 4. This result is robust to the statistical variations that occur when drawing a different instance of the same turbulent density distribution. -In comparison to dynamical friction in a homogeneous medium, presented in Ostriker (1999), we report that for a decelerating BH in a turbulent medium, magnitude of the force depends only weakly on M as long as M > 1. Overall, the r max required to fit the formula from Ostriker (1999) is unrealistically high for high M, where we find that r max exceeds the size of the box probed here by several orders of magnitude. By contrast, r max < ∆x for low M at low values of turbulent wavenumber k, but approaches physically reasonable values of a few pc as k → ∞. We therefore report that in a turbulent medium, dynamical friction is more efficient for high Mach numbers, but less efficient as the BH approaches the transonic regime (i.e. as M → 1). -At a given Mach number, scatter in dynamical friction magnitude is larger for turbulence with lower k than for higher k, due to the more extended but less frequent density perturbations traversed by the BH. There are a few caveats to the work presented here. One is that we do not capture the exchange of momentum between the BH and the background gas due to accretion. Using the BHL accretion ratė we estimate that the mass accretion rate during the supersonic regime with v rel = 10 km s −1 , v rms = 4 km s −1 is ∼5 M Myr −1 . This leads to about 10 3 M accreted during the 200 Myr of supersonic regime, and, hence, a decrease of approximately a tenth of the initial BH velocity due to gas accretion. This effect could potentially allow the BH to slow down into the sub-sonic regime, as accretion becomes more efficient at low Mach numbers. Including the impact of accretion in the long-term evolution of the BH is left to future work. On the other hand, the aftermath of gas accretion, involving injection of energy as feedback, has been shown to act counter to dynamical friction. Indeed, with radiative (Park & Bogdanović 2017), thermal (Souza Lima et al. 2017), and kinetic (Gruzinov et al. 2020) feedback the wake can be destroyed, nullifying the deceleration. Furthermore, feedback can create an underdense region behind the BH and the net effect is for the BH to experience an acceleration in the direction of motion, at least temporarily and for some gas densities and BH masses (Toyouchi et al. 2020). We caution that the small size of our box means that the BH traverses the same region of gas repeatedly throughout the study, despite choosing an angle that minimises this effect. The turbulent forcing implemented here re-randomises the gas density field before the next pass, but we cannot rule out that cumulative effects might affect the late-time evolution of our BH.
Overall, we conclude that dynamical friction due to a turbulent background medium is able to efficiently slow a BH down to M = 1, but cannot further decelerate the BH. This force is self-consistently captured by hydrodynamical simulations as long as the minimum resolution criterion of r BH /∆x > 5 is met, but under-resolved otherwise. When resolved, the force in a turbulent medium exceeds that in a homogeneous medium at high Mach number, but is reduced around the trans-sonic regime.

Appendix B: Dynamical friction in a smooth medium
Creating a smooth equivalent to our turbulent boxes is not trivial for several reasons. Firstly, without the reheating from turbulence, the gas in the box would bulk-cool over the timescales required for the dynamic study. When the gas is not allowed to cool, it is instead dynamically heated by the repeated passage of the BH, and also picks up a non-negligible bulk velocity through momentum transfer from the BH. Both effects artificially speed up the decay of the Mach number, invalidating measurements of the drag force as a function of Mach number. One solution is to use a bigger simulation box, as this increases the thermal and dynamical mass of the gas, but maintaining constant resolution for bigger boxes quickly becomes numerically unfeasible. Through experiment, we found that doubling the box size of our fiducial simulation, Run0, allowed for a reasonable compromise: for an 80 pc box, a simulation with the same input parameters as Run0 (including resolution) sees a final sound speed increased by a factor ∼1.5 and a final bulk speed of the gas of only 0.6 km/s. There is also the problem that the gas does not stay smooth over the timescales required for our setup: old wakes from the repeated passage of the BH through the box create density fluctuations encountered by the BH during its next passage. To tackle these issue, we use a set of simulations that have the same parameters as Run0, but have initially uniform density (and no turbulent forcing), an 80 pc box and a starting Mach numbers of Dimensionless drag force f on the black hole for a suite of non-turbulent equivalents of Run0 with different initial Mach numbers, and the points of our k = 16 simulation for comparison. Black lines represent the analytic force according to Ostriker (1999) for the same values of Coulomb logarithm than in the Fig. 11 for the k = 16 case: 0.1, 3.6, 8. M ini = 12.5, 10, 7.5, 5, and 2.5 respectively. We evolve each simulation only for ∆M = 3, that is for a limited change in Mach number.
Comparing the dimensionless drag force on the BH for our non-turbulent boxes (Fig. B.1) to the highest wavenumber from the turbulent set (bottom right panel of Fig. 11), the evolution of the force as a function of the Mach number in a smooth medium is similar to highest wave-number case, k = 16, but it has a higher magnitude at low M. This is consistent with the previous trend of higher forces at low Mach number for increasing k, and suggests that the suppression of the drag force in the transonic regime is due to the turbulence. The same is true for r max shown in Fig. B.2, which show large values at high Mach numbers and reasonable results of a few pc in the trans-sonic regime. r max = 5 would mean ln(λ) = 3.44, in reasonable agreement with both Chapon et al. (2011;ln(Λ) = 3.2) and Beckmann et al. (2018; ln(Λ) = 3.9). There is little existing work that numerically tests the analytic prediction for the magnitude of the drag force for a uniform medium at high Mach numbers, so it is difficult to tell if the fact that the drag force at high M is up to 50% higher than in the predicted analytic case for constant Coulomb logarithm is a generic result or not. Beckmann et al. (2017) do have a simulation at M = 10 and results possibly also show a higher drag force at such high Mach numbers than the analytic prediction, but do not explore this further. This might be an interesting phenomenon to be explored in future work.
Firmer conclusions are difficult to draw due to the large scatter, induced by the BH crossing its own wake during repeated passages of the box.