Morphology and Mach Number Dependence of Subsonic Bondi–Hoyle Accretion

We carry out three-dimensional computations of the accretion rate onto an object (of size R sink and mass m) as it moves through a uniform medium at a subsonic speed v ∞. The object is treated as a fully absorbing boundary (e.g., a black hole). In contrast to early conjectures, we show that for an accretor with Rsink≪RA=2Gm/v∞2 in a gaseous medium with adiabatic index γ = 5/3, the accretion rate is independent of Mach number and is determined only by m and the gas entropy. Our numerical simulations are conducted using two different numerical schemes via the Athena++ and Arepo hydrodynamics solvers, which reach nearly identical steady-state solutions. We find that pressure gradients generated by the isentropic compression of the flow near the accretor are sufficient to suspend much of the surrounding gas in a near-hydrostatic equilibrium, just as predicted from the spherical Bondi–Hoyle calculation. Indeed, the accretion rates for steady flow match the Bondi–Hoyle rate, and are indicative of isentropic flow for subsonic motion where no shocks occur. We also find that the accretion drag may be predicted using the Safronov number, Θ = R A /R sink, and is much less than the dynamical friction for sufficiently small accretors (R sink ≪ R A ).


INTRODUCTION
Over many decades, the problem of Bondi-Hoyle-Lyttleton (BHL) accretion flow has been studied through analytical estimates and numerical calculations in both two and three dimensions, beginning with Hunt (1971).These studies have covered a range of adiabatic indices and have treated stationary accretors as well as those with either subsonic or supersonic velocities, aiming to determine the rates of mass and momentum accretion as well as the shock structure (if present).
BHL flows are relevant to many astrophysical problems, particularly to the accretion onto black holes.Other accreting objects have hydrodynamically relevant surfaces which limit their accretion -such as protoplanets within disks -though BHL may still be used to obtain an upper bound for the accretion rate (Nelson & Benz 2003).Wind accretion onto stars and X-ray binaries is a commonly-used application of BHL flows, though this typically involves supersonic motion and is therefore not relevant to this work (see, for example, Struck et al. 2004).On the other hand, galaxies in clusters can move through the intergalactic medium (IGM) subsonically due to the high temperature of the IGM (Edgar 2004).Thus subsonic BHL flow may describe accretion onto galaxies from the IGM, as observed by De Young et al. (1980), though galactic outflows may complicate this process (Owen et al. 2000).Finally, binary interactions such as common envelope evolution involve the subsonic or supersonic motion of an accretor through a stellar envelope.Here the axisymmetry of BHL flow is broken by density gradients in the envelope, creating an angular momentum barrier to accretion.Thus the accretion rate must be modified according to the density scale height of the envelope (Taam & Bodenheimer 1989;MacLeod & Ramirez-Ruiz 2015).These applications vary in terms of their velocities, accretion radii, and adiabatic indices (γ), necessitating a general formulation for the accretion rate.
Several attempts have been made to determine a self-similar solution for the accretion flow and to derive an analytical form for the accretion rate Ṁ .However, they have struggled to provide generality in γ and to bridge the gap between stationary, subsonic, and supersonic accretors.The solution for spherical flow of Bondi (1952) is well-known, as is the accretion rate for zero-temperature axisymmetric flow (Hoyle & Lyttleton 1939), which is relevant to highly-supersonic velocities.For intermediate Mach numbers 0 < M ≲ 2, the Bondi rate is often interpolated such that it reproduces these results in the limits M → 0 and M → ∞, respectively (Bondi 1952).Additionally, Ruffert (1994) obtained fitting formulae from the results of their numerical simulations, and subsequently Foglizzo & Ruffert (1997) derived an analytical interpolation between subsonic and supersonic flows for γ > 9/7.These predictions agree with their numerical results for M = 0.6 and 1.4, and suggest that the accretion rate is equal to the spherical Bondi rate for at least some range of Mach numbers.However, the behavior of Ṁ over the full range of subsonic flow is not known.
The relevant length scales are the accretion radius, defined in Hoyle & Lyttleton (1939) as as well as the Bondi radius Here m is the mass of the accretor, v ∞ is its velocity relative to the ambient medium, and c s,∞ is the sound speed.We primarily use R A throughout this work for consistency between our simulations and for comparison to previous works (Ruffert & Arnett 1994;Ruffert 1994;MacLeod & Ramirez-Ruiz 2015).The astronomical objects listed above vary in the ratio of their physical radii R to R A , ranging from R/R A ∼ 1 for galaxies and planets down to 10 −5 -10 −10 for black holes.As there is no a priori reason to believe that the accretion flow is described by the same solution in both regimes, numerical calculations are needed to determine which accretion flows may be described by the BHL model.In this paper, we aim to characterize the flow around both large and small accretors relative to R A , and in particular to determine the accretion rate.We organize this paper as follows.We review the relevant theory to determine the accretion rate in Section 2, and describe our numerical setup in Section 3. We analyze our results in Sections 4 and 5 including the accretion rate, drag forces and location of the sonic point.We discuss these results and conclude in Section 6.

Spherical Accretion
Though we consider axisymmetric flow in this work, the streamlines deflect toward the radial direction due to gravity near m.As the streamlines converge around the accretor, transverse pressure gradients enforce radial flow.We first assume that the flow is one-dimensional in the radial direction and that it has reached a steady state.For subsonic velocities, there are no shocks and the flow is isentropic.
We now review the theory of spherical hydrodynamic accretion and direct interested readers to chapter 14 of Shapiro & Teukolsky (1983) and chapter 13 of Zeldovich & Novikov (1971) for in-depth analysis.Consider a fluid with density ρ, fluid velocity v, Mach number M, and pressure P flowing through a pipe with cross-sectional area A. In the classical theory of nozzles and diffusers (with body forces absent), a relation between the cross-sectional area and fluid velocity can be obtained from conservation of mass and momentum: This area-velocity relation demonstrates that the sign of the acceleration of the flow depends both on the sign of dA and whether the flow is subsonic or supersonic.If a gravitational potential is included in the conservation of momentum, a modified version of the area-velocity relation is obtained: Here the sign of dv depends not only on dA but on dA relative to dΦ, meaning that the flow can potentially be accelerated to supersonic velocities upstream of the throat due to gravity.Spherical accretion onto a sink region can be viewed in the same light as a nozzle in which the area is given simply by A = 4πr 2 and the throat has area A sink = 4πr 2 sink , with a potential Φ = −Gm/r.Then dA = 8πrdr and dΦ = Gmdr/r 2 , so Equation (4) becomes (5) To find the sonic point, we set M = 1, yielding We determine the sound speed at the sonic line, c sonic , from the Bernoulli constant In steady flow, B is constant along each streamline (Chapman 2000), which allows us to equate B sonic = B ∞ .Using the Bernoulli constant in conjunction with (6), we arrive at This is similar to the calculation of Shapiro & Teukolsky (1983) with the caveat that our B ∞ contains a kinetic term v 2 ∞ /2 as in Foglizzo & Ruffert (1997).For γ = 5/3, this gives r sonic = 0, and for γ > 5/3 it gives r sonic < 0. Only for γ < 5/3 is r sonic positive.However, this calculation neglected the finite radius of the accretor.Within R sink , the pressure and density are negligible, which sets a minimum radius for the sonic surface at r = R sink .
For steady radial flow, the mass accretion rate can be expressed as Ṁ = Aρv = 4πr 2 ρv.Because the flow is isentropic, P = Kρ γ for some pseudoentropy K, and c 2 s = γP/ρ = γKρ γ−1 .Combining this with (6) gives the density at the sonic line Placing this and the sound speed into Ṁ , we obtain Ṁ = 4π γK Gm 2γK For γ = 5/3, the prefactor simplifies to However, we encounter the issue that r sonic = 0 for γ = 5/3.This is alleviated by the fact that Finally, this leads us to We see that for the special case of γ = 5/3, Ṁ is independent of r sonic and of the Mach number.This rate is equal to that found by Bondi (1952) for spherically-symmetric accretion from a gas at rest at infinity, which we refer to hereafter as ṀB .

Axisymmetric Accretion Rates
In addition to the accretion rate assuming locally spherial flow discussed above, there are several theoretical predictions and numerical fits for Ṁ .These include: (i) The oft-quoted formulation proposed by Bondi (1952) as an interpolation between the accretion rates for static and supersonic accretors.
The original formula in Bondi (1952) was smaller than ( 14) by a factor of 2, but this was corrected in subsequent works so that in the limit v ∞ ≫ c s,∞ it agreed with the Hoyle & Lyttleton (1939) Shima et al. 1985;Ruffert & Arnett 1994).
(ii) A set of fitting formulae proposed by Ruffert (1994) based on their numerical results (their Equations 2 -8).These formulae were fit to simulations performed at M = 0.6, 1.4 and 10.Interestingly, though only one Mach number below unity was studied, this fit predicts (like Equation 13) that Ṁ is independent of Mach number for subsonic accretors.
(iii) The analytical interpolation formulae for axisymmetric flows with γ > 9/7 proposed by Foglizzo & Ruffert (1997) (their Equations 107 and 108): which describes Ṁ in terms of deviations from the spherical Bondi rate ṀB .Here the entropy is characterized by the "effective Mach number" M e , defined as with the caveat that M e ≥ 1, where 0.5 < λ < 1 is a parameter which can be tuned such that ṀFR97 approaches the Hoyle-Lyttleton rate for high M ∞ .If the flow is taken to be isentropic (M e = 1), then for γ = 5/3 we recover ṀFR97 = ṀB .
In this work, we aim to compare these predictions with the results of numerical calculations.To that end, we now discuss the hydrodynamics solvers with which we tackle this problem.

NUMERICAL SETUP
We perform calculations using the arbitrary Lagrangian-Eulerian code Arepo (Springel 2010) and the Eulerian code Athena++ Stone et al. (2020).A comparison of these numerical solvers and mesh structures is summarized in Table 1.Though they are both grid codes, they differ markedly in their handling of the mesh.Arepo uses a moving, unstructured mesh, utilizing a Voronoi tessellation (Okabe et al. 1992) to draw a grid around an arbitrary distribution of points.These mesh-generating points then move at the local fluid velocity, with corrections to ensure roundness of the cells.Additionally, we impose refinement conditions for the maximum cell size as a function of distance from the accretor.We initialize the Arepo grid with 10 sub-grids, each with 100 cells, where the smallest subgrid's size is determined such that there are at least 1000 cells surrounding the accretor.Our Athena++ calculations use a static spherical-polar mesh in which the cell spacing in the polar and azimuthal directions is constant, with n θ = 40 and n ϕ = 10.In the radial direction, we use n r = 800 cells whose radial extent is geometric, with each cell differing in radial size by a factor of 1.01 from its neighbors.
Both codes handle time integration via a second-order van Leer predictor-corrector scheme, and spatial reconstruction with a second-order piecewise-linear method (PLM).We set the Courant-Friedrichs-Lewy number to η CFL = 0.3 or 0.4.Slope limiters are applied using the van Leer limiter with the corrections described in Mignone (2014) to keep the limiter total-variation diminishing.Additionally, both codes use the Harten-Lax-van Leer contact (HLLC) approximate Riemann solver (Toro et al. 1994).We have tested the effects of other Riemann solvers and slope limiters such as the Roe solver (Roe 1981) in Athena++, the exact iterative solver in Arepo and the Superbee limiter (Roe 1986) in Arepo.We also tested the effect of gas self-gravity in Arepo with R sink /R A = 0.02 and 0.002 but found no significant change to the accretion rate nor the flow morphology in any of these cases.This is due to our choice of physical parameters leading to a gas mass inside R A to always be a few orders of magnitude smaller than the Jeans mass, as well as the accretor.Therefore, in the scenarios checked here, the gravity is always dominated only by the accretor.We note that the self-gravity of the gas is expected to be important in cases of high density regime, as in common envelope evolution.

Outer Boundary
Our wind-tunnel setup requires gas of a specified state to inflow from one end of the domain and freely exit the other end.In Athena++, in which we use a spherical domain, this is accomplished by separately treating the upstream (θ < π/2) and downstream (θ > π/2) sections of the outer boundary.The upstream section imposes the specified ambient conditions in ghost cells, while the downstream portion is a diode: gas is allowed to outflow with zero-gradient conditions but is not allowed to enter the domain.
Arepo uses a cubical domain, handling inflow via an injection region on the upstream end of the box.Mesh cells are created with the specified fluid state within this injection region of width 3 in our internal units (defined in Section 3.2), where at each time-step all cells within this region are updated for the properties of the ambient flow.Outflow is handled by deleting cells whose mesh-generating points are about to leave the domain.By default, the length of the simulation domain is chosen as D = 32R A to match Ruffert & Arnett (1994).In Arepo, this refers to the side length of the cubical domain, while in Athena++ it refers to the diameter of the spherical domain.

Inner Boundary
In both codes, the inner accretion boundary is imposed by creating a region of radius R sink in which the density and pressure are much lower than their surroundings.In Athena++, the inner boundary of the domain is at R sink .
Table 2. Input parameters and derived quantities for all simulations performed in this work: the accretor radius R sink , the type of inner boundary condition, the specific heat ratio γ, the incoming Mach number M∞, the gravitational parameter Gm, the accretion radius RA, the width of the simulation domain D, the mass accretion rate Ṁ , the predicted accretion rate ṀB from (13) (or (10) for γ = 4/3), the ratio of Ṁ to its lower limit Ṁgeo, and the momentum accretion rate Facc.All quantities are expressed in code units, as described in Section 3.2.Immediately inside of this boundary is a shell of ghost cells with state (ρ, v, P ) = (10 −2 ρ ∞ , ⃗ 0, 10 −2 P ∞ ).We take a different approach in Arepo, instead imposing a "sink" region within R sink with a sink particle at its center.Any cell whose vertex falls within this region is eliminated, and its previous mass and momenta are added to those of the sink particle.

Physical Parameters
Unless otherwise specified, we model the gas as ideal with specific heat ratio γ = 5/3.Following Ruffert & Arnett (1994), we choose our code units such that ρ ∞ = c s,∞ = 1, which necessitates P ∞ = K ∞ = 1/γ.Ruffert & Arnett (1994) also chose the characteristic length scale as R A = 2Gm/v 2 ∞ = 1, which sets the (initial) mass of the accretor at Gm = v 2 ∞ /2.We deviate from this somewhat by varying v ∞ while holding Gm constant at the fiducial value of (0.6) 2 /2 = 0.18 to study the dependence of Ṁ only on v ∞ .This means that it is not always possible for the simulation domain to span the entire accretion radius since R A → ∞ as v ∞ → 0. However, the domain size always far exceeds the Bondi radius R B = 2Gm/c 2 s,∞ = 0.36.For completeness, we also perform several runs at at different value of the accretor mass Gm = 0.045.Additionally, we have checked that our results are converged by increasing the linear resolution of our fiducial simulation by a factor of 4 and confirming that our principal results -the accretion rate and accretion drag -are not altered significantly.Ruffert (1994) found that the rates of accretion of both mass Ṁ and momentum F acc are dependent on the size of the accretor R sink , and their numerical results suggest that Ṁ converges for sufficiently small R sink (see their Fig.11).Indeed, the assumptions of Section 2 hold provided there exists some radius within which the flow is radial, spherically symmetric, and steady.Motivated by this, we also vary the accretor size from R sink /R A = 1 down to 0.0002, serving both to confirm the convergence of Ṁ and to determine the value to which it converges.To our knowledge, this is the smallest accretor which has been tested, with Ruffert (1994) using a minimum accretor size of R sink /R A = 0.02 and MacLeod & Ramirez-Ruiz (2015) using R sink /R A = 0.01.
Although the analytics of Section 2 are applicable to the flow near an accretor with R sink ≪ R A and γ = 5/3, it is also relevant to consider larger accretors and other specific heat ratios.For this reason, we also conduct suites of simulations with R sink = R A .With accretors of this size, it is important to properly model the flow in the far field.Thus, for these runs we vary the domain size such that D = 32R A is preserved.We also perform two runs with γ = 4/3, in which the sonic line is predicted by Equation ( 8) to occur at a finite radius.
A summary of all of our simulation parameters can be found in Table 2.The fork of Athena++ used to perform the simulations in the paper is located at https://github.com/ljprust/athena/tree/windtunnel,configured with the problem generator file src/pgen/windtunnel.cpp.The parameter file can be found at inputs/hydro/athinput.windtunnel sink and can be modified according to Table 2 to reproduce our results.

FLOW MORPHOLOGY
Our simulations eventually settle into a steady state, as defined by the variations in derived quantities.Slice plots of the Mach number in both codes for R s = 0.02R A , M = 0.6 are shown in Fig. 1 for γ = 5/3 (left panels) and γ = 4/3 (right panels).The overall morphology is very similar between the two codes.The flow accelerates as it nears the accretor, achieving M = 1 at R sink for γ = 5/3 and at some r > R sink for γ = 4/3.These features are discussed below, as are the accretion rates and drag forces.

Hydrostatic Equilibrium
The gas is compressed as it approaches the accretor, creating pressure gradients which may be sufficient to enforce a degree of hydrostatic equilibrium.Conservation of momentum can be expressed as so the condition for hydrostatic equilibrium (∇P/ρ ≈ −∇Φ) is that the kinetic term is relatively small.We compute each term in Equation ( 17) along the radial direction and plot them along with their sum in Fig. 2 for R sink /R A = 0.002, M = 0.6.We see that ∇P/ρ ≈ −∇Φ for r ≲ R A with the exception of small radii r ≲ 10 −2 R A , where pressure support is lost and the flow accelerates.This indicates that much of the flow near the accretor is suspended in near-hydrostatic equilibrium, where pressure gradients in the radial direction are sufficient to prevent gas from free-falling into the accretor.Gradients transverse to the radial direction are also present and may enforce one-dimensional radial flow, as we now discuss.

Reaching Spherical Symmetry
We can learn if the flow near the accretor is radial and spherically-symmetric by checking the flow properties at the surface of the accretor.In Fig. 3, we show these properties as well as the Mach number transverse to the radial direction M ⊥ and mass accretion rate per unit area Ṁ /dA for the layer of cells adjacent to the sink.We see that for R sink = 0.002R A , M ∞ = 0.9 (left panel) all properties are uniform, indicating spherical symmetry has been reached.The exception is the transverse Mach number, though it is small compared to the radial Mach number M = 1 and decreases with accretor size.A large accretor with R sink = R A , M ∞ = 0.9 exhibits a different morphology; spherical symmetry is not obeyed, and both the radial and transverse flow is supersonic over a portion of the surface.We find violation of spherical symmetry for M ∞ ≳ 0.6 only in the case of R sink = R A , whereas symmetry is obeyed for R sink ≤ 0.02R A regardless of Mach number.
We compare the radial profiles of several quantities for R sink /R A = 0.002 in Fig. 4 along the upstream (left panel) and downstream (right panel) directions in both Arepo and Athena++.The profiles are very similar between the two codes with the exception of the region very near to the accretor surface, which is likely due to the different methods of implementing the inner boundary.The stagnation point is evident in the downstream Mach number profiles at r ≈ 0.07R A .The Mach number reaches unity near R sink from both sides, as predicted in Section 2, while the entropy remains constant at its initial value K ∞ = 1/γ.

Sonic Line
As discussed in Section 2, as gas flows toward the accretor there must exist a point at which it reaches sonic velocity, which for γ = 5/3 occurs at the surface of the accretor r sonic = R sink .We find this to be the case in our simulations, as shown in Fig. 4. For γ < 5/3, the sonic line occurs at a finite radius outside of R sink as estimated from (8).Specifically, for γ = 4/3 we obtain r sonic = 0.0425R A , which is overplotted in the lower right panel of Fig. 1.We see that the predicted and actual sonic surfaces roughly agree, though the actual sonic surface is not spherically symmetric.Because the region within the sonic surface is causally disconnected from the outside, the accretion rate cannot be affected by the flow within the sonic surface.This means that, unlike γ = 5/3, the accretion rate is unaffected by the size of the accretor provided R sink ≲ r sonic .Indeed, for γ = 4/3 we find no dependence of Ṁ on R sink (see Table 2).Here the only uniform quantity is entropy, with supersonic flow transverse to the accretor over a portion of its surface.All quantities are presented in code units, as defined in Section 3.2.

Stagnation Point
Downstream of the accretor, there exists a neutral point in the flow where M = 0.The position of this stagnation point cannot be determined analytically, but an upper bound for its distance R stag from the accretor can be estimated.Lyttleton (1972) argued that R stag may be of order the accretion radius, as it separates streamlines which will escape

Accretion Drag
The material passing through the accretor surface carries not just mass, but momentum as well.When the flow is nearly radial, much of this momentum cancels out, though small asymmetries can lead to a net force.In Arepo, the momentum accretion is computed by taking the momentum of all cells which enter the sink region and adding to that of the sink particle.As the sink particle is allowed to move in Arepo, it acquires a small but nonzero velocity ≈ 2 × 10 −8 in the lab frame for R sink /R A = 0.002 as a result of accreting an amount of gas equal to ≈ 10 −5 of its mass.In Athena++, the static mesh prevents such motion, though the Arepo results demonstrate the validity of this approximation.Here we compute the rate of momentum accretion along the symmetry axis by again integrating over the surface of the accretor: where the results are shown in Table 2 and in Fig. 7. F acc is generally positive (denoting a force in the direction of the freestream flow), with some exceptions.This can be compared with the gravitational drag due to dynamical friction, which for a subsonic accretor takes on the form 1999).For our chosen parameters ρ ∞ = c s,∞ = 1, Gm = 0.18, and M ∞ = 0.6, the dynamical friction is F DF ≈ 0.1.Thus, we can infer from Fig. 7 that dynamical friction is the dominant source of drag for small accretors with R sink ≲ 10 −1 R A .
Though the accretion drag depends on the accretor size R sink , unlike the mass accretion rate we do not find convergence for a sufficiently small sink region.For a non-gravitating body, one would expect the cross-section for accretion drag to be σ = πR 2 sink .On the other hand, Safronov (1972) suggests that for a gravitating accretor the cross-section is given by σ = πR 2 sink (1 + Θ), where Θ = R A /R sink is the Safronov number.For Θ ≪ 1, this produces the expected behavior σ = πR 2 sink , while in the limit Θ ≫ 1 it reduces to σ = πR 2 A Θ −1 .Here, we find that the resulting drag well-approximates the data, as shown in Fig. 7.This is valid only for γ = 5/3, as for other equations of state the cross-section is altered due to the sonic surface.

DISCUSSION AND CONCLUSIONS
In this paper, we have studied Bondi-Hoyle-Lyttleton accretion using the Athena++ Eulerian code as well as the moving-mesh code Arepo.We assume axisymmetry and isentropic flow, which is relevant to subsonic accretors.These codes give nearly identical results despite their many differences in both the numerical evolution of the fluid and in the measurement of derived quantities, indicating that our results are robust to the numerical scheme used.We show that the gas within the accretion radius is suspended in hydrostatic equilibrium, with the exception of a narrow shell adjacent to the accretor, where the gas accelerates to sonic velocity.For γ = 5/3, this results in an accretion rate which is independent of Mach number and is equal to the spherical Bondi rate (Bondi 1952), depending only on the accretor mass and the gas pseudoentropy.This result suggests a physical explanation in which the flow near to the accretor surface can be approximated as radial and spherically-symmetric.For subsonic flow, there are no shocks and therefore no change in entropy.Previously this Mach number invariance has not been explicity confirmed, resulting in many works still using the interpolation of the Bondi rate Ṁ ∝ (v 2 ∞ + c 2 s,∞ ) 3/2 for subsonic flow.We contrast this result with several previous predictions.We find that the interpolated Bondi rate ( 14) (Bondi 1952) overestimates Ṁ for M < 1, while our rate is consistent with the interpolation of Foglizzo & Ruffert (1997) only for M ≲ 0.7.Our results are consistent with the fitting formulae of Ruffert (1994) -which is also devoid of Mach number dependence in the subsonic regime -while offering a much simpler calculation.
For γ > 5/3, the accretion rate ( 13) is not valid due to the vanishing of the sonic radius.However, for γ < 5/3 the sonic surface occurs at a finite radius outside of the accretor, given by Equation (8).Indeed, the case γ = 4/3 is of interest when the gaseous medium is radiation-pressure dominated.Our runs with γ = 4/3 are consistent with theory: the actual location of the sonic surface agrees reasonably well with (8), and the Ṁ corresponding to this sonic radius agrees with the rate measured from our simulations to within several percent.Additionally, unlike the case of γ = 5/3, the accretion rate is independent of accretor size provided the accretor is enclosed by the sonic surface.
This rate is applicable to accretors without hydrodynamically relevant surfaces (e.g.black holes) moving through fairly homogeneous gaseous media.This occurs, for example, when black holes accrete from the interstellar medium, or in binary stellar processes such as common envelope evolution.In particular, if the accretor is engulfed by a star, both the convective and radiative regions of stellar structure may be treated by the formalism in Section 2. Fortunately, the γ values that cannot make use of this formalism, γ > 5/3, are rarely relevant for accretion flows -though see Richards et al. (2021) for a discussion of Bondi accretion with a stiff equation of state in the context of neutron stars.
We find a net force exerted on the accretor due to the acquisition of momentum from accreted gas.For γ = 5/3, this accretion drag is in the direction of the freestream flow and its value may be predicted using the Safronov number.However, analytical estimates show that the force from accretion drag is eclipsed by that of the dynamical friction when R sink ≪ R A .
Simulations of larger accretors with R sink = R A exhibit several qualitative differences from small ones.Here Ṁ far exceeds the spherical Bondi rate and approaches the lower bound Ṁgeo , which is the limit in which gravity is ignored and matter is simply swept up by the motion of the accretor.For sufficiently high Mach number, the radial, spherically-symmetric solution is absent, instead exhibiting supersonic flow in both the radial and transverse directions.This indicates that accretion onto objects with radii comparable to their accretion radius cannot be modeled by the formalism of Section 2. In this paper, we have treated the gas as homogeneous, neglecting gradients which are expected to be present in many astrophysical scenarios.These gradients could create an angular momentum barrier to accretion or violate our assumption of radial flow near the accretor.Turbulent media are not explored here, and may affect the results.We have also not considered supersonic accretors, which would result in the accretion of high-entropy shocked material.According to (13), this would reduce Ṁ .We leave the study of supersonic accretors to future work.

Figure 2 .Figure 3 .
Figure 2. Comparison of the terms in the momentum equation (17) normalized by the gravitational acceleration for R sink = 0.002RA, M∞ = 0.6 in Athena++.For intermediate radii 10 −2 ≲ r/RA ≲ 1, the pressure gradients are comparable to the potential, indicating pseudo-hydrostatic equilibrium.

Table 1 .
Comparison of the numerical algorithms in our simulation codes.Numerical values reported are from our R sink = 0.002RA runs.