Gas Dynamical Friction on Accreting Objects

The drag force experienced by astronomical objects moving through gaseous media (gas dynamical friction) plays a crucial role in their orbital evolution. Ostriker derived a formula for gas dynamical friction by linear analysis, and its validity has been confirmed through subsequent numerical simulations. However, the effect of gas accretion onto the objects on the dynamical friction is yet to be understood. In this study, we investigate the Mach number dependence of dynamical friction considering gas accretion through three-dimensional nested-grid simulations. We find that the net frictional force, determined by the sum of the gravitational force exerted by surrounding gas and momentum flux transferred by accreting gas, is independent of the resolution of simulations. Only the gas outside the Bondi–Hoyle–Lyttleton radius contributes to dynamical friction, because the gas inside this radius is eventually absorbed by the central object and returns the momentum obtained through the gravitational interaction with it. In the subsonic case, the front–back asymmetry induced by gas accretion leads to larger dynamical friction than predicted by the linear theory. Conversely, in the slightly supersonic case with a Mach number between 1 and 1.5, the nonlinear effect leads to a modification of the density distribution in a way that reduces the dynamical friction, compared with the linear theory. At a higher Mach number, the modification becomes insignificant and the dynamical friction can be estimated with the linear theory. We also provide a fitting formula for dynamical friction based on our simulations, which can be used in a variety of applications.


INTRODUCTION
Moving objects interacting with their surrounding media, such as stars and gases, are ubiquitous in the universe.Such moving objects include binary stars orbiting in common envelopes (Ricker & Taam 2008), planets orbiting in protoplanetary disks (Teyssandier et al. 2013), and supermassive black holes (SMBHs) with > 10 6 M ⊙ moving in galaxies and in some cases forming pairs with other SMBHs (Narayan 2000;Armitage & Natarajan 2005;Escala et al. 2004Escala et al. , 2005)).
The formation of SMBH binaries resulting from galaxy mergers has been of great interest for a long time.After two galaxies hosting their own SMBHs merge, the SMBHs are believed to reduce their distance through Corresponding author: Tomoya Suzuguchi suzuguchi@tap.scphys.kyoto-u.ac.jp interactions with stars and gases until they eventually merge emitting gravitational waves (GWs) (Begelman et al. 1980).Recently, some SMBH binaries separated on the kpc scale have been observed as dual active galactic nuclei (see, e.g., De Rosa et al. 2019, for a review).Moreover, NANOGrav has presented compelling evidence of a stochastic GW background consistent with SMBH binary mergers (Agazie et al. 2023).
Recent discoveries of SMBHs at redshifts greater than six have prompted investigations into their origins in the early universe (see, e.g., Inayoshi et al. 2020, for a review).While various SMBH formation scenarios have been proposed, most of them presume that SMBHs have grown from smaller BHs, predicting a large number of intermediate-mass black holes (IMBHs) with ≲ 10 6 M ⊙ harbored in primeval galaxies.Consequently, mergers of such galaxies lead to formation of IMBH binaries, whose GWs emitted at their mergers are a major target for space interferometers such as Laser Interferometer Space Antenna (LISA, Amaro-Seoane et al. 2017, 2023), TianQin (Luo et al. 2016), and Taiji (Ruan et al. 2020).
To estimate the event rate of the GWs from IMBH binaries, understanding their cosmological evolution is crucial.However, it is often impractical to follow the evolution of IMBH binaries in cosmological simulations resolving pc or sub-pc scale dynamics, where they interact with surrounding stars and gases.Therefore, it is necessary to introduce an appropriate sub-grid model that can describe the effect of such interactions on the motion of IMBH binaries in cosmological simulations (e.g.Beckmann et al. 2019).The interactions of a moving object with its surrounding medium have been studied for a long time.Chandrasekhar (1943) explored the concept of dynamical friction, which is a resistive force that acts on an object moving in a stellar distribution.He developed an analytical formula to describe this force by considering the interactions between the object and the stars as a two-body gravitational problem.Later, Ostriker (1999) (hereafter O99) studied the gas dynamical friction, which is the resistive force resulting from the movement of an object through a surrounding gaseous medium (see also Just & Kegel 1990).She examined the linear density perturbations caused by the object and derived a formula for the frictional force by integrating the gravitational force from the surrounding gas.
The analytical formula proposed by O99 has been extensively tested through numerical simulations, and it has been confirmed that the formula accurately predicts the resistive force in both linear and circular motion cases, as long as the moving object only causes small density fluctuations, for which linear analysis is applicable (e.g.Sánchez-Salcedo & Brandenburg 1999, 2001;Kim & Kim 2007;Escala et al. 2004;Kim et al. 2008).However, when a density perturbation of order unity appears in the vicinity of the object, the formula obtained by the linear analysis deviates from the actual friction force in the nonlinear regime (e.g.Kim & Kim 2009;Thun et al. 2016).It is noteworthy that these previous simulations did not take into account the accretion of gas.Further numerical research is needed to gain a better understanding of the effects of gas accretion on the gas dynamical friction.
Gas accretion onto moving objects has been studied through analytic consideration (e.g.Hoyle & Lyttleton 1939;Bondi & Hoyle 1944) and hydrodynamic simulations (e.g., Hunt 1971Hunt , 1979;;Shima et al. 1985;Ruffert 1994a,b;Ruffert & Arnett 1994;Ruffert 1995Ruffert , 1996;;Edgar 2004 for a review) for the case of uniform linear motion in a homogeneous medium.The case in which the central object is a BH and exerts radiative feedback on the surrounding gas has also been studied through radiative hydrodynamic simulations (e.g. Park & Ricotti 2013;Park & Bogdanović 2017;Toyouchi et al. 2020;Sugimura & Ricotti 2020;Ogata et al. 2021Ogata et al. , 2024)).However, the focus of the above studies is mainly on mass accretion growth rather than the gas dynamical friction1 .Understanding of the latter is still limited even in the simplest setup of linear motion in homogeneous medium without any feedback.While Shima et al. (1985) and Ruffert (1995Ruffert ( , 1996) ) have obtained the resistive force acting on the objects, taking into account the effects of gas accretion, the force apparently depends on the resolution, and the numerical convergence of dynamical friction has not been addressed.Furthermore, their simulations were limited to a few specific Mach numbers, and the Mach number dependence of dynamical friction has not been fully determined.Therefore, a comprehensive comparison with Ostriker's formula is yet to be done in cases with gas accretion.
In this paper, we perform three-dimensional numerical simulations to investigate the dynamical friction exerted on an object moving at a constant velocity in a homogeneous gaseous medium while accreting the surrounding gas.We calculate the dynamical friction as a net force due to the gravitational pull of surrounding gas and momentum flux of accreting gas.We find that the net force converges with increasing spatial resolution, as the resolution dependences of the gravitational pull and momentum flux cancel out with each other.Furthermore, through simulations with varying Mach numbers, we obtain the dynamical friction as a function of the Mach number, including the full nonlinear effects with accretion, and give a physical interpretation by comparing our results with Ostriker's linear theory.We provide a dynamical friction formula based on our results, which is an updated version of Ostriker's formula and can serve as a more realistic subgrid model in future cosmological simulations in which the gas dynamical friction cannot be directly calculated.
The rest of the paper is organized as follows.Section 2 provides a brief overview of the theoretical basis, including the linear perturbation theory of the gas dynamical friction given by O99.Section 3 outlines our numerical simulations and the cases examined.Section 4 presents the results of the numerical simulations and offers a physical interpretation of the results.Section 5 describes the application of our friction formula and compares our results with those of previous studies.Finally, Section 6 summarizes the paper.

THEORETICAL BACKGROUND
Here, we provide a brief overview of the theoretical background that is useful in interpreting our simulation results.

Characteristic length scales
Consider a situation in which a static object accretes the surrounding gas in a homogeneous medium.We assume the polytropic equation of state with an index γ for the gas.In this situation, we can obtain a solution for transonic accretion flow with spherical symmetry (Bondi 1952).The characteristic length scale of the flow, known as the Bondi radius, is given by where M and c s is mass of the object and sound speed, respectively.At a radius sufficiently smaller than R B , the accretion flow can be approximated with free fall; the gas density and inward velocity follow ρ ∝ r −3/2 and u ∝ r −2 , respectively.The accretion rate is given by where ρ ∞ is the gas density at infinity, and λ(γ) is a factor of order unity that depends on γ (e.g.Shapiro & Teukolsky 1983).For example, in the isothermal case (γ = 1), λ = e 3/2 /4.Next, consider a case in which the object has a relative velocity with respect to the surrounding homogeneous medium, and the accretion flow becomes axisymmetric.If the relative velocity is highly supersonic (v ≫ c s ), the gas with the impact parameter smaller than the Hoyle-Lyttleton radius (Hoyle & Lyttleton 1939;Bondi & Hoyle 1944), is accreted onto the object.When the relative velocity is comparable to the sound speed (v ∼ c s ), the typical length scale is often expressed with the Bondi-Hoyle-Lyttleton (BHL) radius, which is an interpolation of the Bondi and Hoyle-Lyttleton radii.

Gas dynamical friction in the linear perturbation theory
In this subsection, we briefly review the gas dynamical friction formula derived by O99.The formula is based on linear analysis and assumes a gravitational source with mass M moving at a constant velocity v that appears at time t = 0 in a homogeneous medium with density ρ ∞ .The gas is assumed to follow the isothermal equation of state with c s .Under these assumptions, the time evolution of linear density perturbations excited by the moving object is analytically solved.The integration of the gravitational force due to the density perturbation yields the dynamical friction, F DF , as a function of the Mach number M ≡ v/c s .
The resulting formula is written as with a M-dependent factor the Coulomb logarithm and a normalization factor with the size of radial momentum flux in the spherical (Bondi) accretion The radius r min denotes the minimum radius of the gas that contributes to the dynamical friction.In the supersonic case (M > 1), r min was introduced as an artificial lower limit for the radial integration of the gravitational force to avoid the divergence in the integration (O99) and cannot be determined within the framework of linear theory.Later in this paper, we determine r min based on our simulation results, such that it reproduces our numerical results according to Equation (6).Note that the above formula depends on t when M > 1, but not when M < 1.In the linear theory, the frictional force is roughly proportional to ln (r max /r min ), where r max is the maximum radius of gas contributing to the dynamical friction.In the subsonic case, the radial extent to which the density perturbation propagates gives r max ≃ (c s + v)t, while the radius within which the density distribution maintains a front-back symmetry gives r min ≃ (c s − v)t.Therefore, the frictional force becomes time independent with the time dependence in ln (r max /r min ) canceling out.However, in the supersonic case, r max is again ≃ (c s + v)t but r min does not have a time dependence.Therefore, in contrast to the subsonic case, the time dependence in ln (r max /r min ) does not cancel out, and the dynamical friction becomes time dependent.
Finally, we make three remarks regarding the linear formula given by Equation (5).Firstly, the formula is not valid near M = 1, as Equation ( 5) diverges logarithmically at M = 1.Furthermore, in the supersonic case, the formula is derived under the assumption that M > 1+r min /(c s t) (O99).Secondly, when gas accretion is present, the linear approximation is accurate only in regions far away from the object (r ≫ R BHL ) where the density perturbation is small.Lastly, the formula does not take into account the effect of accretion.

Simulation setup
We consider a body moving linearly in a uniform isothermal gaseous medium.We use the Cartesian coordinates (x, y, z), with −x being the direction of motion, and the velocity of the body can be written as v = −ve x , where e x denotes the unit vector in the x direction.In our simulations, we adopt the frame of reference of the body, i.e., we set a uniform initial velocity u(x, t = 0) = −v to the gas with the position of the body fixed, where u(x, t) is the gas velocity at the position x and the time t.Once a simulation starts, the gravitational force of the body begins to pull surrounding gas towards the body.The body is treated as a sink particle that accretes gas within the radius R sink when the density exceeds the critical value n cr .The gravitational softening is effective only within the sink radius (Matsumoto et al. 2015).
Our fiducial parameter set corresponds to the case of an IMBH moving through a dense medium, a situation possibly realized in the early universe (e.g., Milosavljević et al. 2009b,a;Park & Ricotti 2012, 2013;Inayoshi et al. 2016;Sugimura et al. 2017;Inayoshi et al. 2022;Sugimura et al. 2018;Takeo et al. 2018;Toyouchi et al. 2019Toyouchi et al. , 2020Toyouchi et al. , 2021; see also Inayoshi et al. 2020 for a review).However, our results are applicable to a wide variety of cases as long as the radius of objects is sufficiently smaller than Bondi radius (see Equation 1) because the hydrodynamic equations we solve are essentially scale-free.We set the initial gas density to n 0 = n ∞ = 10 5 cm −3 and the critical density for the sink accretion to n cr = (1/50)n 0 .As long as the simulations resolve the supersonic free-falling part of the accretion flows, i.e., R sink ≪ R BHL , the value of n cr does not affect the flow structures outside R sink .2Considering the warm neutral medium, we set the gas temperature to T 0 = 10 4 K and the mean molecular weight to µ ≃ 1.26.These values give the gas sound speed The mass of the object is set to M = 10 3 M ⊙ , and the corresponding Bondi radius is Note that the increase in the body's mass due to accretion during a typical simulation period of t = 10 6 yr (see Section 3.2) is approximately ∆M ∼ ṀB t ∼ 2 × 10 3 M ⊙ and comparable to the body's mass.However, we neglect the mass increment of the body in this work, in order to develop an understanding of gas dynamical friction in a steady state.

Cases considered
We explore a variety of cases by changing the Mach number M of the moving IMBH, including both subsonic and supersonic cases.The details of these cases are summarized in Table 1.We continue a simulation run until the accretion rate becomes nearly constant.The final time t final is listed in Table 1.Note-Col.1:Mach number of the relative velocity, Col.2: ratio of the Bondi-Holye-Littleton radius to the Bondi radius, Col.3: size of the computational domain, Col.4: radius of the sink particle, Col.5: highest AMR refinement level, Col.6: time by which the accretion rate converges to a nearly constant value.In all cases, we set the mass of the object to M = 10 3 M⊙, the initial gas density to n0 = 10 5 cm −3 , and the gas temperature to T0 = 10 4 K.
With the AMR module in the code, we set nested grids in the computational box, placing higher resolution grids near the central object.The number of the AMR base grids is set to (N x , N y , N z ) = (128,128,128) for all the cases.We take the larger size of the computational domain R box (defined as half the side length of computational box) with the larger Mach number for the subsonic cases.This is because, when the Mach number is higher, it takes more time for the accretion rate to reach a steady state, which is determined by R B /|v − c s | for M < 1.The sink radius is chosen to sufficiently resolve the BHL radius in each calculation.Since the BHL radius becomes smaller as the Mach number increases, we set R sink ≃ 1.95 × 10 3 au for M < 2.0, and take values that are half and a quarter of this value for M = 2.23 and 2.73, respectively.The minimum grid size in each calculation, R grid , is set to be one-fourth of the sink radius.The highest level of AMR,

Evaluation of dynamical friction
We evaluate the dynamical friction exerted on the moving body in the following manner.From a simulation snapshot, we calculate the gravitational pull of surrounding gas in the direction of the body's motion as where V is the entire computational domain excluding the interior of the sink particle.Next, we compute the momentum flux transferred by accreting gas in the di-rection of the object's motions where S sink is the surface of the sink particle.We then obtain the net frictional force, which we simply call dynamical friction throughout this paper, as the sum of these quantities, According to the above definition, the body is decelerated when F DF < 0 and accelerated when F DF > 0.

4.1.
Representative cases with M ≃ 0.6 and 1.5 Here we focus specifically on the cases with M ≃ 0.6 and 1.5 as representative subsonic and supersonic cases, respectively.We analyze the time evolution of the accretion rate and friction, paying attention to numerical convergence with improving spatial resolution.

Snapshots of density perturbations
Figure 1 presents the temporal evolution of density perturbations induced by a body moving at a subsonic speed of M ≃ 0.62.Panels (A) to (C) describe the time evolution, illustrating the outward expansion of the region influenced by the central body.It is evident that the downstream side of the body exhibits higher density, leading to a gravitational force acting in the direction of deceleration.Figure 2 is the same as Figure 1 but for a supersonic case with M ≃ 1.49.The influence region has a conical and spherical shape due to the supersonic speed of the body (see Figure 2 of O99).Similarly to the subsonic case, the downstream side of the body has a higher density, and the gravitational force acts in the direction of deceleration.
Figure 3 shows the same snapshots as in Figures 1 (upper) and 2 (lower), but for enlarged views of the accretion flow near the body.This figure shows that gas is pulled by the gravitational force of the accretor and moves towards it within the BHL radius, both in subsonic and supersonic cases.Outside the BHL radius, the gas motion is still affected by the accretor's gravity, but the gas passes by the accretor without being trapped by it.

Time evolution of accretion rate and friction
The top panel of Figure 4 displays the time evolution of the accretion rate.In both subsonic and supersonic cases, the accretion rates take nearly constant values after approximately 10 5 yr, indicating that the accretion flows reach quasi-steady states.
The bottom panel of Figure 4 describes the time evolution of the frictional force.In this panel, we plot the net frictional force F DF (Equation 13), together with the gravitational force by surrounding gas F grav (Equa- tion 11) and the momentum flux transferred by accreting gas F acc (Equation 12).We put minus signs on F DF and F grav , as F DF and F grav are negative (causing deceleration) while F acc is positive (causing acceleration).
In the subsonic case, F grav approaches a constant value, which is consistent with the linear perturbation theory (see Section 2).The momentum flux F acc also becomes constant after the accretion rate becomes con-  2), while the vertical axis in the lower panel represents the forces exerted on the body normalized by FB = 4π(GM ) 2 ρ∞/c 2 s (see Equation 8).The solid line represents the net frictional force −FDF, which is the difference between the gravity force from density perturbations −Fgrav (dashed) and the momentum flux of the accreting gas Facc (dot-dashed line).Note that −FDF becomes constant for t ≳ 10 5 years in the subsonic case, and it continues to increase in the supersonic case.
stant.As a result, the net frictional force, F DF = F grav + F acc , converges to a constant value in the quasisteady state.Hereafter, we consider this converged value as dynamical friction in the subsonic case.
In the supersonic case, −F grav increases with time in proportion to log t, as predicted by the linear theory (see Section 2).On the other hand, F acc becomes constant after the gas accretion rate becomes constant, as in the subsonic case.Consequently, −F DF increases with time, reflecting the time dependence of F grav .8).The different colors represent different values of M, where the blue and red lines correspond to M ≃ 0.62 and ≃ 1.49, respectively.We evaluate the gravitational and net frictional forces after they reach constant values in the case with M ≃ 0.62 and at t = 1.0 × 10 5 yr in the case with M ≃ 1.49.The square symbols represent our fiducial choice of the sink radius, R sink = 0.14RB (Table 1).The vertical arrows indicate the positions of the BHL radii with M ≃ 0.62 and 1.49.For each case of M ≃ 0.62 and 1.49, while −Fgrav(Facc) increases (decreses) monotonically with decreasing sink radius, −FDF converges to a constant value for R sink ≪ RBHL.
To assess the numerical convergence of the frictional force with increasing resolution, we examine the impact of varying the size of the sink particle on the magnitude of friction.Figure 5 shows the magnitudes of the gravitational force and the dynamical friction as functions of the sink radius.We see that the magnitude of gravitational force increases as the sink radius decreases, whereas the magnitude of dynamical friction remains nearly constant.
The gravitational force F grav is affected by the size of the sink.According to Section 2, the gas inside the BHL radius is approximately in free fall, and thus the density profile can be expressed as ρ ∝ r −3/2 , where r is the distance from the body.Ignoring the angular dependence of the radial distribution of the density, we can roughly estimate |F grav | as This means that, as the sink size decreases, the body experiences a stronger gravitational force from the highdensity gas closer to it.As a result, F grav is resolution dependent and cannot be considered a physical force.The cause of net frictional force being independent of the sink radius is as follows.The gas inside the BHL radius is mainly absorbed by the body, transferring momentum to it.This momentum counteracts the gravitational force from the gas that is enclosed in the BHL radius.Consequently, the sum of the gravitational force and the momentum flux does not depend on the sink radius.We consider this combined force to be the dynamical friction.We also verified that the dynamical friction converges to a constant value for other cases with M ≃ 1.29, 1.73, 2.23, and 2.73, although this is not shown in the figure for simplicity.

Mach number dependence of dynamical friction
In what follows, we investigate the Mach number dependence of the dynamical friction, considering all the cases with different Mach numbers (Table 1).
Figure 6 illustrates the variation of the net frictional force (dynamical friction) with different Mach numbers.In the range of M > 1, different colors of the data points represent different epochs of 63 t B (5 × 10 5 yr, orange), 125 t B (1 × 10 6 yr, blue), and 250 t B (2 × 10 6 yr, green), where For the subsonic cases, except for M ≃ 0.99, the black dots represent the values in the quasi-steady state.We also compare our simulations with the linear perturbation theory (O99) with the time fixed to t = 125 t B .In Figure 6, the blue dashed line represents the linear perturbation theory, for which we choose r min = R BHL /2 to give an overall fit to the simulation data (see Equation 5-8). Figure 6 shows that on the subsonic side, the simulation results tend to be higher than those of linear perturbation theory.On the supersonic side, the simulation results are smaller than those predicted by the linear perturbation theory.However, the difference decreases as the Mach number increases.
We provide a formula that fits our simulation results much better than the linear perturbation theory with r min = R BHL /2 but is based on it,  8) and the horizontal axis represents the Mach number M. The different colors denote different periods: orange, blue, and green represent t = 63 tB (5 × 10 5 yr), 125 tB (1 × 10 6 yr), and 250 tB (2 × 10 6 yr), respectively, with black corresponding to time-independent quantities.Here, tB = RB/cs is the sound crossing time at the Bondi radius (see Equation 15).In subsonic cases, except for M ≃ 0.99, we plot the dynamical friction at 25 tB (2 × 10 5 yr) for M ≃ 0.19, 0.37, 0.62 and at 63 tB (5 × 10 5 yr) for M ≃ 0.87, when the dynamical friction has already reached a constant value.We overplot our fitting formula given by Equation ( 16) for t = 63, 125, and 250 tB (solid lines) and the linear formula from O99 given by Equation ( 5) with rmin = RBHL/2 for 125 tB (dashed line).
where I(M) takes distinct forms in different regimes of M. First, in the subsonic regime of 0 < M < 0.99, where The functional shape of Equation( 17) is the same as that of linear theory Equation (6) for M < 1 except for the factor of λ(M).We introduce this factor to enhance the frictional force with respect to the linear values.Our simulation results show the largest deviation from the linear theory around M = 0.5, where λ(M) takes its maximum value.Next, in the supersonic regime of M > 1.1, where This formula also has the same form as that of Equation ( 6 where and The linear function of M in the intermediate regime (Equation 21) connects the function in the subsonic (Equation 17) and supersonic (Equation 19) regimes, making the fitting formula provided by Equation ( 16) a continuous function of Mach number across the entire Mach number range.
The above results suggest that, in the supersonic regime, we can fix the undetermined parameter in the linear theory r min , by applying the linear formula to reproduce the simulation results.Equations ( 19) and ( 20) indicate that r min determined in such a way should be similar to the BHL radius, a comparison that will be given in Section 4.2.2.Here, we consider how our results are compared with Sánchez-Salcedo & Brandenburg (1999) and Kim & Kim (2009), who also conducted numerical simulations to study the dynamical friction exerted on the linearly moving body.Note that both of these authors assumed a gravitational softening radius R soft significantly larger than the Bondi radius, neglecting the effect of mass accretion, unlike ours.As a result, they confirmed that their results are well described by the linear theory assuming r min sufficiently larger than the Bondi radius.For example, Kim & Kim (2009) gave r min = 0.35M 0.6 R soft , which corresponds to 35 R B ≤ r min ≤ 68 R B within the range of 1.0 ≤ M ≤ 3.0.This is in stark contrast to our results, where the gravitational softening radius is much smaller than the Bondi radius, and the influence of mass accretion is of central interest.
So far, we have presented the overall Mach number dependence of dynamical friction.In the subsequent sections, to understand the obtained Mach number dependence, we examine how the dynamical friction is caused in our simulations comparing our results with the linear perturbation theory.We describe the subsonic case in Section 4.2.1 and the supersonic case in Section 4.2.2.

Subsonic cases
First, we consider why the value of dynamical friction in the subsonic case is larger than that obtained from the linear perturbation theory (see Figure 6).In Figure 7 (I), the isodensity contours in the linear perturbation theory form symmetric elliptical shapes centered at the body.In Figure 8, we plot the gas number density profiles along the x-axis in both forward and backward directions with respect to the body's motion.The profiles show that the symmetric structure observed in Figure 7 (I) extends up to x = (c s − v)t on the x-axis.Consequently, in the linear theory, the gravitational force exerted by the gas for x < (c s − v)t cancels out due to backward-forward symmetry and friction is caused by the asymmetry in the gas density for x > (c s − v)t.
In contrast, in Figure 7 (II), the centers of the isodensity contours in our simulation are located downstream side of the body.This positional shift occurs because accretion of gas onto a body moving at a constant velocity takes place primarily on the downstream side, leading to an accumulation of gas in the downstream direction.Thus, the density on the downstream side be-  ) shows the density predicted by the linear perturbation theory (O99), while Panel (II) shows that from our simulation.The white filled circle at the center illustrates the sink particle and white lines indicate the isodensity contours between 2×10 5 cm −3 and 7×10 5 cm −3 with intervals of 10 5 cm −3 (solid) and between 1.8 × 10 5 cm −3 and 1.4 × 10 5 cm −3 with intervals of 2 × 10 4 cm −3 (dashed).The plotted region has a side length of 6 × 10 4 au, and Panel (II) corresponds to a three-times zoom-in of panel (C) in Figure1.The orange and blue dashed lines represent the positions of the Bondi and BHL radii, respectively.comes higher than on the upstream side.In Figure 8, the density profiles along the x-axis exhibit asymmetric even in the region with x < (c s − v)t, where the gas distribution is symmetric in the linear theory.Hence, in our simulation, the gravitational force from the region of x < (c s − v)t does not cancel out, unlike in the linear theory as mentioned above.This explains why the dynamical friction is greater than that predicted by the linear theory.Recall that the gas within x < R BHL does not contribute to dynamical friction because mass accretion compensates for the gravitational force.Only the gas in the range R BHL < x < (c s − v)t provides an 10 3 10 4 10 5 10 6 x [au] This work Ostriker (1999) Figure 8. Density profiles along the x-axis for the case with M ≃ 0.62 at the time of Figure1 (C).The horizontal axis represents the distance x from the central body, while the vertical axis represents the normalized difference between the density at x and its ambient value, i.e., (n − n∞)/n∞.We show the density obtained from our simulations with blue lines and that obtained from the linear perturbation theory (O99) with orange lines.Solid and dashed lines represent values for x > 0 and x < 0, respectively.Vertical dotted lines indicate, from left to right, the BHL radius RBHL, (cs − v)t, and (cs + v)t.The region shaded in gray represents the location of a sink particle.Note that in the linear perturbation theory, the density distribution is symmetrical for |x| < (cs − v)t, leading to identical density profiles for x > 0 and x < 0.
additional contribution to dynamical friction, compared to the linear theory.

Supersonic cases
We have used the identical functional form from the linear theory in our fitting formula for the supersonic cases (Equation 19), albeit with the substitution of r min with ζ(M)R BHL .We note that, in the linear theory, r min represents the minimum radius of the gas that contributes to the dynamical friction.Using the linear theory formula to fit our nonlinear simulation results, the factor ζ(M)R BHL accounts for additional effects beyond what is considered in the linear theory.Here we extract such additional effects by considering the Mach number dependence of ζ(M).
To do so, we introduce the critical radius r crit , the radius only defined with a given simulation snapshot, apart from the linear theory.This is the minimum radius of the gas that contributes to the dynamical friction.In other words, the gravitational force from the gas enclosed by r crit cancels out with the momentum flux carried by the accreting gas.We calculate r crit from a given simulation snapshot as follows.We compute the gravitational force in the direction of the body's motion Figure 9.The Mach number dependence of characteristic length scales in our simulations for the supersonic cases.
The blue circles represent the critical radius rcrit, where the momentum flux carried by the accreting gas cancels out the gravitational force exerted by the gas within a radial distance r, FDF(< r) = 0 (Equation 26).The orange squares represent the minimum radius r min,fit , which is determined so that the linear formula (Equation 6) fits the dynamical friction measured in our simulations, FDF (Equation 13).These radii divided by RBHL are plotted as functions of the Mach number M in the figure.The orange curve illustrates ζ(M) = rmin/RBHL given by Equation ( 20), which is a fitting function for the Mach number dependence of r min,fit .The horizontal dotted line represents r/RBHL = 1/2.
exerted by the gas within a radial distance r where V (r) represents the spherical volume of the radius r.We consider the dynamical friction only from the gas within the volume V (r) and obtain r crit as the radius at which Figure 9 shows how the radius r crit varies as a function of the Mach number M. We see that r crit is always comparable to R BHL for the entire range examined of 1 < M < 3. The gas within the BHL radius R BHL accretes in free fall and returns its momentum to the body, which does not contribute to the friction.This coincides with the definition of r crit , which explains why r crit ∼ R BHL .
By comparing ζ(M) and r crit /R BHL , we can analyze the differences between our simulation results and the predictions of linear theory, which correspond to the  The left-hand side of Figure 10 presents the density dis-tributions obtained from the linear perturbation theory (upper) and our simulation (lower) at M = 1.06.In the linear perturbation theory, the density distribution takes the form of a conical plus spherical shape (see Figure 2 of O99).In contrast, the density distribution in the simulation exhibits a different shape, with the surface of density discontinuity located on the upstream side of the body.In this case, the presence of high-density regions on the upstream side results in a smaller density difference between the upstream and downstream sides.Consequently, the dynamical friction is reduced compared to the linear theory, leading to a larger value of ζ(M) compared to r crit /R BHL (Figure 9).The flow structure at M = 1.06 resembles that of the subsonic case (see Figure 1).We propose that the M = 1.06 case lies in the transitional regime between the subsonic and supersonic cases.Consequently, the friction resulting from the gas outside r crit is comparable in both the linear theory and our simulation, leading to ζ(M) ∼ r crit /R BHL (Figure 9).We observe that the flow structure for M ≥ 2 exhibits more temporal variability compared to other cases due to the instability reported in the literature (e.g., Matsuda et al. 1987;Fryxell & Taam 1988).Figure 11 shows the temporal changes in the gas distribution for the case of M = 2.23.These temporal variations result in fluctuations in r crit .However, after t ≥ 1 × 10 5 yr, r crit stabilizes at a constant value, with fluctuations of approximately 10%.To obtain the values of r crit shown in Figure 9, we average the results over six different time intervals, spaced roughly every 1 × 10 4 yr within the range of (1.0 − 1.5) × 10 5 yr.

Application of friction formula
Our fitting formula for the dynamical friction, Equation ( 16), can be applied to a variety of problems.It is an updated version of the formula proposed by O99, which has been widely used in the literature.One example of its application is in modeling the motion of a BH within a galaxy, as observed in cosmological simulations with limited spatial resolution.For instance, Dubois et al. (2013) employ the formula where α (≥ 1) represents a boost factor that accounts for the unresolved dense gas in the vicinity of the BH, ρ and cs denote the mean density and mean sound speed of the resolved gas, and I lin (M) is the same function as defined in Equation (6).Our formula can also be applied to the above equation by substituting I lin (M) with I(M) (Equation 16).One of its primary effects is to systematically enhance the friction force within the subsonic range, compared to the linear theory presented in Equation ( 6).
In the supersonic range, the function I(M) is time-dependent.Dubois et al. (2013) set ln Λ = ln(r max /r min ) = 3.2 based on the simulation results given by Chapon et al. (2013).Since our simulations have shown that r min should be taken as ζ(M)R BHL ,3 one can determine r max by some criteria.A simple choice of r max is to take the typical size of the system, such as that of a galactic disc.When considering a body in linear motion, as in our simulations, one can substitute r max = Mc s t.For a body in circular motions, such as in a binary system, r max can be approximated as the orbital radius (e.g.Kim & Kim 2007).However, since Kim & Kim (2007) only focus on the linear regime, the dynamical friction of accreting bodies in circular motions still needs to be addressed in future work.Shima et al. (1985) conducted a seminal numerical study on the friction experienced by an accreting body in linear motion at a constant velocity.Building upon this, Ruffert (1996) performed a more comprehensive analysis using 3D hydrodynamic simulations with a nested grid technique.The latter work considered various cases where a body moving at different velocities (M = 0.6, 1.4, 3.0, and 10) accretes a nearly isothermal gas (with an adiabatic index of γ = 1.01).However, the main focus of this study was to determine the accretion rates onto the body rather than to evaluate the frictonal  Ruffert (1996) and those obtained from our formula (Equation 16), respectively.These symbols correspond to the results for the cases of M = 0.6 at 54 tB (4.3 × 10 5 yr; orange), M = 1.4 at 29 tB (2.3 × 10 5 yr; green), and M = 3.0 at 2 tB (1.6 × 10 4 yr; red).The purple dashed line exhibits the formula provided in Escala et al. (2004), which does not have explicit time dependence.

Accreting body in linear motion
force.We here compare Ruffert (1996) with our findings.We have carefully examined the time dependence and numerical convergence with improving the spatial resolution of the net frictional force.
Figure 12 illustrates the Mach number dependence of dynamical friction, similar to Figure 6.However, in this case, we compare our formula (Equation 16; represented by stars) with the results obtained by Ruffert (1996) (squares).Note that Ruffert (1996) followed the evolution for different periods depending on the Mach numbers: 54 t B (4.3×10 5 yr) for M = 0.6, 29 t B (2.3×10 5 yr) for M = 1.4,and 2 t B (1.6 × 10 4 yr) for M = 3.0.To obtain our corresponding results, we substitute these times into our formula for the supersonic cases with M = 1.4 and 3.0.There is no time dependence for the subsonic case with M = 0.6.Our results are consistent with the overall trend observed in Ruffert (1996).However, there is a discrepancy of approximately 20% between the data points at M = 1.4.This discrepancy can likely be attributed to the limited size of the computational domain used in Ruffert (1996), which does not fully encompass the entire region affected by density perturbations.As a result, the friction, which should have increased steadily over time, remains constant in their study, leading to an underestimation of dynamical friction.

Non-accreting body in a circular orbit
We also compare our results with previous simulations on the dynamical friction of a body moving in a circular orbit without accretion.Specifically, we compare our results with Escala et al. (2004), who performed SPH simulations of circularly orbiting binary BHs within a non-singular isothermal sphere.They derived a formula to describe the velocity dependence of dynamical friction based on their simulation results.
In Figure 12, we show the formula provided by Escala et al. (2004) with the purple dashed line.In the subsonic regime, the friction calculated with the formula of Escala et al. (2004) appears to be underestimated compared with our results, which might be attributed to their exclusion of gas accretion in the simulations.In the supersonic case, their formula lacks explicit time dependence, unlike ours, possibly because the region contributing to friction is determined by the size of the core of the non-singular isothermal sphere in their simulation, as discussed in Section 5.1.However, a direct comparison is not straightforward due to other numerous differences between their and our simulations, including the treatment of accretion (included or ignored) and the assumed trajectories (linear or circular).We point out that their formula may be specialized to specific situations, which is complementary to our formula designed for general use.
For non-accreting bodies, Sánchez-Salcedo & Brandenburg ( 2001) and Escala et al. (2004) argue that, around M ≃ 1, the dynamical friction for circular orbits is suppressed compared to that for linear motion given by Sánchez-Salcedo & Brandenburg (1999), who reproduce the results of O99.In future studies, we intend to investigate the dynamical friction of accreting bodies that are in circular orbits.

CONCLUSIONS
We have investigated the dynamical friction acting on accreting bodies moving at a constant velocity through 3D nested-grid simulations, incorporating gas accretion with sink particles.We have systematically explored the velocity dependence of dynamical friction over a broad range of Mach numbers, 0 < M < 3.0.Furthermore, we have compared our results with the linear theory (Ostriker 1999), to examine the impact of gas accretion on the dynamical friction.Our results are summarized as follows: • The net force acting on the body is the sum of gravitational force exerted by the surrounding gas and momentum flux transferred by the accreting gas (we refer to this net force as dynamical friction).Through our simulations, we find that the dynamical friction is independent of the sink radius, while the gravitational force of the surrounding gas increases as the sink radius decreases.The gas between the sink and BHL radii does not contribute to the dynamical friction, because the gas within the BHL radius accelerated by the body's gravity subsequently falls onto the body, returning the momentum gained before falling.Due to the cancellation of the gravitational force within the BHL radius and the momentum flux, the dynamical friction is effectively determined by gravitational force of the gas outside the BHL radius.
• In the subsonic case, the dynamical friction in our simulation is larger than that of linear theory.According to the linear theory, the region just outside the BHL radius has front-back symmetry and does not contribute to the dynamical friction.However, in our simulations, the effect of accretion breaks the front-back symmetry and introduces a new contribution from the gas in this region, leading to larger dynamical friction.
• In the supersonic case, the linear theory O99 encompasses a parameter, r min , representing the minimum radius of the gas contributing to the dynamical friction.Although r min cannot be determined within the framework of linear theory, our simulations of accreting bodies reveal that it is physically determined.Assuming the same functional form as the linear formula (O99), we find that the dynamical friction in our simulations can be reproduced with r min ∼ R BHL due to the cancellation of the gravitational force and the momentum flux within the BHL radius, as described above.
• We also find that r min /R BHL is a decreasing function for 1 < M < 3.0, which converges to become constant at r min /R BHL = 1/2 for M ≥ 3 (Equation 20).We attribute higher r min /R BHL for 1 ≤ M ≤ 1.3, i.e., smaller dynamical friction compared with the linear formula with r min = R BHL /2, to the nonlinear effect that reduces the asymmetry in the density distribution (see Figure 10).In the slightly supersonic case with 1 ≲ M ≲ 1.3, considered as transitional regime from the subsonic to supersonic regimes, the gas distribution has a similarity to the subsonic case, resulting in the largest deviation from the linear theory that predicts a distinct transition between the two regimes.
• We provide the fitting formula of dynamical friction for accreting bodies (Equation 16-23), based on our simulation results.Our formula serves as an updated and improved version of the original linear formula proposed by O99.Similar to the linear formula, our revised formula can be used in a wide range of applications, such as galaxy formation simulations, where it can be employed as a subgrid model for a moving and accreting BH.
In conclusion, we determine the velocity dependence of dynamical friction experienced by astronomical objects moving in a uniform medium, taking gas accretion into account.Note that our simulations are based on idealized setup with uniform gas density and that physical effects such as magnetic fields and radiation feedback are not considered in this study.Nonetheless, our results serve as a crucial initial step towards understanding the motion of astronomical objects, such as BHs, as they traverse a medium while accreting gas.
Figure 1.The evolution of density perturbations caused by a body moving at a subsonic speed of M ≃ 0.62.Panels (A), (B), and (C) represent the time sequence, corresponding to epochs at (A) 6×10 4 yr, (B) 1.2×10 5 yr and (C) 2.4×10 5 yr, respectively.The body is assumed to move horizontally to the left.The white filled circle at the center illustrates the sink particle.The solid lines represent isodensity contours, with the closest to the center denoting 5 × 10 5 cm −3 and the others spaced at intervals of 10 5 cm −3 .The dashed lines are also contours for the less dense outer part, with the closest to the center denoting 1.8 × 10 5 cm −3 and the others spaced at intervals of 0.2 × 10 5 cm −3 .The orange and blue dashed circles represent the Bondi and BHL radius, respectively.
Figure 3. Enlarged views of the accretion flows for cases with M ≃ 0.6 (upper) and M ≃ 1.5 (lower).We take the snapshots at 2.4 × 10 5 yr (the same as panel (C) in Figure 1) for the upper panel and 1.0 × 10 5 yr for the lower panel.The red arrows indicate the gas velocity vectors, with the legend in the upper right of each panel showing their normalization.The solid line closest to the center denotes 1 × 10 6 cm −3 and the others spaced at intervals of 2 × 10 5 cm −3 .The dashed lines are also contours for the less dense outer part, with the closest to the center denoting 1.8 × 10 5 cm −3 and the others spaced at intervals of 0.2 × 10 5 cm −3 .The orange and blue dashed circles represent the Bondi and BHL radii, respectively.In the lower panel, we present the result from the run with the sink particle size reduced by a factor of eight compared with the fiducial run (see Section 4.1.3),to provide a better description of the inner flow structure.

Figure 4 .
Figure 4.The time evolution of the accretion rate (upper) and the friction force (lower).In both panels, the blue and red lines represent the same subsonic and supersonic cases as in Figures 1 and 2, with Mach numbers of M ≃ 0.62 and M ≃ 1.49, repsectively.The crosses on these lines correspond to the epochs (A)-(C) presented in these figures.The vertical axis in the upper panel represents the accretion rate normalized by the Bondi rate ṀB (see Equation2), while the vertical axis in the lower panel represents the forces exerted on the body normalized by FB = 4π(GM ) 2 ρ∞/c 2 s (see Equation8).The solid line represents the net frictional force −FDF, which is the difference between the gravity force from density perturbations −Fgrav (dashed) and the momentum flux of the accreting gas Facc (dot-dashed line).Note that −FDF becomes constant for t ≳ 10 5 years in the subsonic case, and it continues to increase in the supersonic case.

Figure 5 .
Figure5.Dependence of the gravitational force from the gas −Fgrav (dashed line), the momentum flux carried by the accretion gas Facc (dot-dashed line), and the net frictional force −FDF = −(Fgrav + Facc) (solid line) on the sink radius.The horizontal axis represents the sink radius normalized by the Bondi radius RB, and the vertical axis represents the force exerted on the body normalized by FB = 4π(GM ) 2 ρ∞/c 2 s (see Equation8).The different colors represent different values of M, where the blue and red lines correspond to M ≃ 0.62 and ≃ 1.49, respectively.We evaluate the gravitational and net frictional forces after they reach constant values in the case with M ≃ 0.62 and at t = 1.0 × 10 5 yr in the case with M ≃ 1.49.The square symbols represent our fiducial choice of the sink radius, R sink = 0.14RB (Table1).The vertical arrows indicate the positions of the BHL radii with M ≃ 0.62 and 1.49.For each case of M ≃ 0.62 and 1.49, while −Fgrav(Facc) increases (decreses) monotonically with decreasing sink radius, −FDF converges to a constant value for R sink ≪ RBHL.

Figure 6 .
Figure 6.Velocity dependence of dynamical friction, obtained by the sum of gravitational force and momentum flux.The vertical axis represents the magnitude of dynamical friction |FDF| normalized by FB = 4π(GM ) 2 ρ∞/c 2 s (see Equation8) and the horizontal axis represents the Mach number M. The different colors denote different periods: orange, blue, and green represent t = 63 tB (5 × 10 5 yr), 125 tB (1 × 10 6 yr), and 250 tB (2 × 10 6 yr), respectively, with black corresponding to time-independent quantities.Here, tB = RB/cs is the sound crossing time at the Bondi radius (see Equation15).In subsonic cases, except for M ≃ 0.99, we plot the dynamical friction at 25 tB (2 × 10 5 yr) for M ≃ 0.19, 0.37, 0.62 and at 63 tB (5 × 10 5 yr) for M ≃ 0.87, when the dynamical friction has already reached a constant value.We overplot our fitting formula given by Equation (16) for t = 63, 125, and 250 tB (solid lines) and the linear formula from O99 given by Equation (5) with rmin = RBHL/2 for 125 tB (dashed line).
) for M > 1 but with the parameter r min replaced by the BHL radius R BHL multiplied by ζ(M) ∼ O(1).The factor ζ(M) decreases with M for 1.1 < M ≤ 3.0 and takes a constant value of 1/2 for the larger Mach numbers.Finally, in the intermediate regime for 0.99 ≤ M ≤ 1.1, Figure 7 compares the gas number density n H near the center in the case of M ≃ 0.62 predicted by the linear perturbation theory (panel (I); O99) and obtained from our simulation (panel (II); at t = 2.4 × 10 5 yr, same as Figure 1 (C)).

Figure 7 .
Figure 7.Comparison of density distributions in the linear theory (panel I) and simulation (panel II) for the case with M ≃ 0.62 at t = 2.4 × 10 5 yr (same as panel (C) in Figure1).Panel (I)  shows the density predicted by the linear perturbation theory (O99), while Panel (II) shows that from our simulation.The white filled circle at the center illustrates the sink particle and white lines indicate the isodensity contours between 2×10 5 cm −3 and 7×10 5 cm −3 with intervals of 10 5 cm −3 (solid) and between 1.8 × 10 5 cm −3 and 1.4 × 10 5 cm −3 with intervals of 2 × 10 4 cm −3 (dashed).The plotted region has a side length of 6 × 10 4 au, and Panel (II) corresponds to a three-times zoom-in of panel (C) in Figure1.The orange and blue dashed lines represent the positions of the Bondi and BHL radii, respectively.

Figure 10 .
Figure10.The same as shown in Figure7, but for the supersonic cases where M ≃ 1.06 (left) and M ≃ 1.49 (right) in an epoch of t ≃ 1.5 × 10 6 yr.Note that for the case with M ≃ 1.49, the result obtained from the run with a sink radius that is half the fiducial size is presented.This run provides a more detailed view of the density structure near the central body, although the value of dynamical friction does not change significantly when the sink radius is halved, as described in Section 4.1.3.nonlinear effects.If ζ(M) is very close to r crit /R BHL , it means that the density distribution outside r crit in our simulation is similar to that predicted by the linear theory.In this case, the nonlinear effects are not significant.This is illustrated in Figure9, where ζ(M) does not completely coincide with r crit /R BHL but the deviation is not significant for 1 < M < 3. The relations r crit ∼ R BHL and ζ(M) ∼ r crit /R BHL explain why ζ(M) ∼ O(1).In particular, for M ≳ 1.3, ζ(M) is very close to r crit /R BHL .For these cases, both ζ(M) and r crit /R BHL are only slightly higher than 1/2.This explains why the dynamical friction obtained in our simulations is slightly lower than the prediction of linear theory when r min = R BHL /2 (see the blue solid and dashed lines in Figure6).Figure 9 illustrates significant discrepancies between ζ(M) and r crit /R BHL for M ≲ 1.3, indicating that the nonlinear effect becomes more prominent as the Mach number decreases below M ∼ 1.3.This behavior of ζ(M) can be explained by comparing the density distributions from our simulations and the linear theory.The left-hand side of Figure 10 presents the density dis-

Figure 9
illustrates significant discrepancies between ζ(M) and r crit /R BHL for M ≲ 1.3, indicating that the nonlinear effect becomes more prominent as the Mach number decreases below M ∼ 1.3.This behavior of ζ(M) can be explained by comparing the density distributions from our simulations and the linear theory.

Figure 11 .
Figure 11.Non-steady flow patterns observed in the case with M ≃ 2.23.We show the time evolution of the density distribution in panels (A)-(C), corresponding to 7.0 × 10 4 yr, 9.0 × 10 4 , yr, and 1.1 × 10 5 yr, respectively.Here, we present the results from the run with R sink = (1/4)R sink,fiducial .The blue dashed circles denote the BHL radius.

Figure 12 .
Figure12.Same as Figure6but for comparison of our results with previous works.The square and star symbols represent the results obtained byRuffert (1996)  and those obtained from our formula (Equation16), respectively.These symbols correspond to the results for the cases of M = 0.6 at 54 tB (4.3 × 10 5 yr; orange), M = 1.4 at 29 tB (2.3 × 10 5 yr; green), and M = 3.0 at 2 tB (1.6 × 10 4 yr; red).The purple dashed line exhibits the formula provided inEscala et al. (2004), which does not have explicit time dependence.