Bondi or not Bondi: the impact of resolution on accretion and drag force modelling for Supermassive Black Holes
Abstract
Whilst in galaxysize simulations, supermassive black holes (SMBH) are entirely handled by subgrid algorithms, computational power now allows the accretion radius of such objects to be resolved in smaller scale simulations. In this paper, we investigate the impact of resolution on two commonly used SMBH subgrid algorithms; the BondiHoyleLyttleton (BHL) formula for accretion onto a point mass, and the related estimate of the drag force exerted onto a point mass by a gaseous medium. We find that when the accretion region around the black hole scales with resolution, and the BHL formula is evaluated using local massaveraged quantities, the accretion algorithm smoothly transitions from the analytic BHL formula (at low resolution) to a supply limited accretion (SLA) scheme (at high resolution). However, when a similar procedure is employed to estimate the drag force it can lead to significant errors in its magnitude, and/or apply this force in the wrong direction in highly resolved simulations. At high Mach numbers and for small accretors, we also find evidence of the advectiveacoustic instability operating in the adiabatic case, and of an instability developing around the wake’s stagnation point in the quasiisothermal case. Moreover, at very high resolution, and Mach numbers above , the flow behind the accretion bow shock becomes entirely dominated by these instabilities. As a result, accretion rates onto the black hole drop by about an order of magnitude in the adiabatic case, compa to the analytic BHL formula.
keywords:
hydrodynamics  accretion  methods: numerical  black hole physics1 Introduction
One common problem in astrophysics is that the range of scales involved in a given problem can be huge, from atomic to galactic, and as such, it is extremely difficult to fully capture a physical process in a single simulation. This is especially true for processes that suffer from an ’inverse cascade’ problem, where (unresolved) small scale behaviour influences the outcome on (resolved) large scales. For these reasons, we have to rely on subgrid models, which aim to replicate the impact of unresolved, small scale behaviour on scales relevant to the simulation at hand, using only information available in the simulation. The success of such a subgrid algorithm therefore depends both on the relevance of the physics it contains, and its ability to account for it over the widest possible range of scales.
Super massive black holes (SMBH) in cosmological or idealised galaxy simulations are one such problem. The gas fuelling the black hole flows from the Mpc scales of the cosmic web, through the kpc scales of the galaxy down to the last stable orbit and event horizon of the black hole on scales. It is currently impossible to adequately track the gas from where it originated all the way to the black hole in a single simulation. Instead, in Adaptive Mesh Refinement (AMR) particlegrid codes such as RAMSES, SMBH are typically modelled as “sink” particles (Krumholz et al., 2004; Dubois et al., 2010), i.e. the black hole is a massive particle that moves over the grid, removing gas from a small accretion region centred on its current location (see Bate et al., 1995; Springel et al., 2005, for the equivalent in Smoothed Particle Hydrodynamics codes). Although there exist alternatives (see next paragraph), accretion usually proceeds at the BondiHoyleLyttleton (BHL) rate. More specifically, the amount of mass accreted by the sink is calculated using the analytic formula of Bondi & Hoyle (1944) and Hoyle & Lyttleton (1939) derived for a point mass accretor moving at constant velocity through an homogeneous background. In this approach, quantities measured “at infinity” should in principle be used to infer the mass accretion rate onto the sink (see Section 2.3 for details).
Whilst a good starting point in situations with limited information, this approach has several notable shortcomings, particularly the question of where “infinity” is to be defined in a busy galaxy simulation. The analytic solution also does not account for the background gravitational potential of the host galaxy (Korol et al., 2016), the selfgravity of the gas, nor does it consider density or velocity gradients, angular momentum or gas effects such as pressure or instabilities (see Edgar (2004) for a detailed review of BHL accretion). A variety of suggestions have been put forward to address some of these shortcomings. They fall broadly into two camps: (i) replace the BHL model with another that takes as input different largescale properties of the host galaxy, such as gravitational torque (Hopkins & Quataert, 2011; AnglésAlcázar et al., 2015), vorticity (Krumholz et al., 2004) or velocity dispersion of the bulge (Hobbs et al., 2011); or (ii) model unresolved physics on smaller scales via accretion disc schemes (Power et al., 2011; Dubois et al., 2014) or unresolved circularisation and viscous transfer considerations (DeBuhr et al., 2010; RosasGuevara et al., 2015). We refer the interested reader to Negri & Volonteri (2017) for a comprehensive review of current accretion schemes.
The BHL analytic solution assumes a uniform density and velocity of the gas at infinity, a situation that most closely resembles simulations where all characteristic lengthscales associated with the SMBH are much smaller than the local resolution, and the gas reservoir per cell is large compared to the accreted mass. Indeed in this situation, “infinity” can reasonably be understood to mean the gas cells surrounding the sink, as the black hole’s gravitational potential is unresolved, and therefore has little impact on the local gas flow. This is particularly true in the absence of feedback, as the cold, dense phase of the interstellar medium (ISM) gas that is expected to feed the black hole, is also underresolved in low resolution simulations. However, while this has led some authors to introduce a simple boost factor to the Bondi rate to compensate for the lack of resolution (Booth & Schaye, 2009), Negri & Volonteri (2017) recently showed that the interplay between local gas density, accretion rate and resolution is much more subtle, and becomes even more so in the presence of black hole feedback.
Computational advances now allow us to resolve a much wider range of scales in galaxy evolution simulations, from the Mpc scales of the cosmic web at large scales to the length scales of the black hole at small scales, fast approaching the physical scale (Schwarschild or Kerr radius) of the black hole itself. In this paper, we investigate how the popular BHL accretion algorithm, where local mass or volume weighted average gas properties are used to approximate quantities at infinity, behaves as the gravitational influence of the black hole becomes better resolved, and the local cell mass becomes comparable to the mass accreted per time resolution element. In that respect, the work presented here is similar to Ruffert (1995a, b), who investigated the validity and convergence of the analytic BHL formula using a range of fixed spherical, inflowonly regions to represent the sink. However, our work focusses less on the validity of the BHL analytic solution itself — even though we do discuss the observed instabilities and how they cause systematic deviations from the analytic solution — than on testing the impact of representing the BH as a sink particle, as well as estimating local gas averages from partially or fully resolved density features within the accretion region of the sink.
Our work also builds on Krumholz et al. (2004), which first introduced the use of Lagrangian sink particles in grid codes, and Dubois et al. (2010), who implemented it in RAMSES. More specifically, we carry out the most thorough exploration of the BHL parameter space performed to date using such a model, both in Mach number and resolution, and for two gas adiabatic indices. This allows us to test the subgrid accretion algorithm behaviour in a variety of regimes were it has or will be used in galaxy simulations, including two features only briefly mentioned in Krumholz et al. (2004): a change in accretion mode at sufficiently high resolution, and accretion becoming inefficient when the resolution approaches the scale radius of the black hole.
As the total mass accreted onto SMBH in galaxy simulations depends not only on the amount of gas removed, but also on the ability of the black hole to stay attached to local high density structures, we then look in detail at the subgrid algorithm used to account for the dynamical friction exerted by the gaseous gravitational wake behind the sink particle. The particular simulations presented here were performed with RAMSES, but we believe that the conclusions we reach are widely applicable to subgrid algorithms where the size of the accretion region, over which local gas properties are measured and from which gas is removed, decreases with increasing resolution.
The structure of the paper is as follows: in Section 2 we present the setup of the simulations, and explain the sink particle algorithm in RAMSES in detail. Section 3 investigates accretion for a range of Mach numbers and resolutions, including a detailed study of the Bondi problem (Section 3.1) and the HoyleLyttleton problem (Section 3.2). Section 4 present results for the impact of Mach number and resolution on dynamical friction, both resolved on the grid and calculated analytically. Section 5 discusses the origin of instabilities observed in the flow. Conclusions can be found in Section 6.
2 The simulations
To study gas dynamics in the vicinity of the sink particle, we set up a series of 3D isolated boxes where the sink particle is embedded in a uniform gas flow. Each simulation is parametrised by a Mach number, , where and are the sound speed and flow velocity respectively. As the sink is held fixed at the centre of the box, is both the absolute flow velocity, and the flow velocity relative to the sink. Resolution is measured as the number of cells, , within the sink particle’s scale radius , so that , where is the size of the smallest grid cell in the simulation.
This characteristic scale radius depends on as follows:
(1) 
where is the gravitational constant and is the mass of the sink particle. Setting for all simulations reduces the number of parameters of the problem. Under this assumption, , so the mass of the sink determines the relative size of the local scale radius to the size of the box. In order to reduce edge effects, we set , which leads to . The final parameter is the gas pressure, which is set to for all simulations. We perform simulations for two values of , an adiabatic one where , and a quasiisothermal one, where .
Most work on BHL accretion in the literature parametrises resolution by , where is the size of the accretor. However, with the sink particle algorithm used here, cannot be unambiguously defined, as we remove gas from a region spanned by a Gaussian kernel, whose width varies from , depending on local conditions (see Section 2.3 for details). We therefore use as the resolution parameter instead, but write i.e. when comparing to previous work.
2.1 Nomenclature
In the remainder of the paper, quantities denoted with a subscript, , are analytical values which parametrise the simulations presented here, defining both initial and boundary conditions. Quantities denoted with a subscript, , are calculated numerically on the fly, and are based on massaveraged quantities within the accretion region of the sink. We explain in detail what this means in Section 2.3. Quantities with a subscript, , are calculated from cell values across the entire simulation box. The term “local” refers to the immediate vicinity of the black hole, whereas the term “global” describes conditions at the boundaries of the simulation, far away from the sink particle. Time averaged values are denoted by triangular brackets, .
The input parameters for each simulation are summarised in its name as follows: a simulation labelled mn is adiabatic () when , and quasiisothermal () when . For example, a simulation called m10.0n4.8a has , and .
Scalar quantities are denoted with italic lettering, , whereas vector quantities are denoted in bold . denotes a radial distance from the accretor, whereas denotes the distance from the accretor measured along the line axis of symmetry of the wake, with positive values measured downstream of the accretor, and negative values upstream.
2.2 Simulation setup
All simulations presented here are performed with the AMR code RAMSES (Teyssier, 2002) and consist of a simple threedimensional isolated cubic box, with a sink particle kept fixed at the centre. Gravity is prescribed by an analytic gravitational field for a point mass, which inherits the mass and location of the sink and a gravitational softening length equal to the smallest cell size, . There is no radiative cooling and the gas is not selfgravitating.
Gas enters the simulation box diagonally, from the bottom left corner, to avoid issues associated with grid aligned flows^{1}^{1}1Such as oddeven decoupling (Quirk, 1994) or the Carbuncle phenomenon at shock fronts (Peery & Imlay, 1988; Elling, 2009). All ghost regions (cells outside of the simulation domain) are set using a zerogradient scheme, except in the case where constant inflow boundary conditions are used instead to avoid edge effect propagation. Accretion proceeds via the usual RAMSES sink particle algorithm (Krumholz et al., 2004; Dubois et al., 2010). We use a linear reconstruction scheme for conservative variables, and a Courant factor of as higher order reconstructions and/or larger Courant factors lead to numerical artefacts.
To reach the required number of levels of refinement within a sufficiently large box, we employ a fixed nested grid strategy, as seen in Figure 1. In all cases, the black hole is surrounded by concentric shells of progressively lower refinement down to a root grid of cells, with the highest resolution shell determined by the chosen resolution value of . Test of various nested grid configurations show that the results presented here are insensitive to the exact grid layout, so the radius of grids was doubled for each level to balance the size of refined regions with computational cost. For all simulations with , all levels of refinement are present from the beginning of the simulation. Simulations with developed a shock feature during the initial settling phase, which led to a permanent breaking of the symmetry of the flow (see Appendix A). This issue can be prevented by running simulations with with a lower initial resolution of until , and refining to the desired level after this preconditioning period. To further avoid numerical effects when refining, only a single additional level of refinement is added per timestep until the correct is reached. In order to give all simulations enough time to damp initial condition transients, all timeaveraged values in the paper are measured for , where the time is given in units of scale radius crossing time:
(2) 
2.3 The accretion algorithm
The standard sink particle accretion algorithm in RAMSES is that implemented by Dubois et al. (2010) from an algorithm originally proposed by Krumholz et al. (2004), using the BondiHoyleLyttleton interpolation formula (Hoyle & Lyttleton, 1939; Bondi & Hoyle, 1944; Edgar, 2004),
(3) 
where the gas properties, such as the density , the bulk velocity and the sound speed are nonlocal measures meant to be taken at infinity, or “far from the gravitational influence of the point mass”.
To calculate these nonlocal properties, RAMSES uses a cloudparticle system (Dubois et al., 2010; Teyssier et al., 2011), where a total of 2109 cloud particles are scattered around the sink with a spacing of , filling a sphere with radius , where is the local cell width. For simulations more complex than the isolated case presented here, the sink is free to move across the grid, making the cloud particles invaluable in sampling local cells, whether the sink is positioned at the cell centre or not. However, in the simulations discussed in this paper, placing the sink at the centre of a cell, combined with the spacing of creates ambiguity in the assignation of parent cell to cloud particles, as some of these particles are located exactly on cell boundaries. RAMSES assigns ambiguous particles preferentially to the cell located downwards and to the left in the coordinate system used here, creating a preferential accretion direction. To avoid this issue, we move all cloud particles radially inward by distributing them with a separation of . This restores spherical symmetry as each cloud particle is unambiguously assigned to a host cell (See Figure 2 for a 2D slice through the cloud of particles).
Each cloud particle samples the properties of its host cell. The massweighted average gas properties used in the BHL formula are then calculated by summing over all cloud particles according to:
(4) 
where the contribution of each cloud particle is weighted by a Gaussian kernel,
(5) 
and is the cloud particle distance from the sink. The width of the kernel, , depends on the value of an interpolated scale radius,
(6) 
where and are the relative velocity and sound speed of the sink host cell respectively. is also limited by the local cell size, so that altogether:
(7) 
Gas mass is removed from local cells, and added to the sink particle, by looping over the cloud particles and removing mass from each host cell according to
(8) 
where is the local timestep. is the accretion rate calculated on the fly according to the BHL formula in Equation 3, evaluated using local massweighted average quantities.
To avoid creating numerical instabilities by removing too much gas mass from a single cell, the total mass accreted per timestep is capped at 75 % of the cloud particle host cell gas mass. This criteria is commonly employed in RAMSES, and ensures that no single cell is emptied in a single timestep when the sink particle moves across the grid. As we force the sink particle to remain in a specific host cell, we also introduce a density floor, , below which accretion is not permitted. We emphasize that this extra parameter is specific to the set of simulations presented in this work rather than a standard RAMSES parameter. Our results are insensitive to the choice of as long as .
2.4 The drag force algorithm
One way in which the black hole interacts with gas in its vicinity is through dynamical friction. The relative velocity between the two creates an overdensity downstream of the black hole, and the enhanced gravitational pull that ensues acts as a drag force. This process transfers momentum from the sink to the gas — contrary to accretion which transfers momentum from the gas to the black hole — thus boosting the tendency of the sink to stay attached to overdense local gaseous structures. When this drag is unresolved, black holes often dynamically decouple from the gas, leading, in the worst cases, to their ejection from the galaxy disc. This has been a recurring problem in galaxy simulations including SMBHs (Sijacki et al., 2007; Volonteri et al., 2016; Biernacki et al., 2017).
RAMSES includes a subgrid model for dynamical friction (Dubois et al., 2013), in which the drag force on the sink is estimated according to the analytic formula derived by Ostriker (1999):
(9) 
where is the unit vector pointing in the direction of the relative velocity between the gas and the sink, and
(10) 
In the previous equation, is the Coulomb logarithm, with the numerical value chosen according to Chapon et al. (2013), is the maximum distance information has travelled downstream of the sink, and is the softening length used to curtail the diverging density profile near the sink. We reexamine the choice for in this paper, as the value from Chapon et al. (2013) is extracted from black hole merger simulations that do not include accretion.
In the simulations presented here, the sink particle is held at a fixed position, so the drag force subgrid algorithm is not employed. However, we compare the analytic formula given by Equation 9 to values measured directly from gravitational wakes of the sink particle when the latter is resolved. As with the BHL accretion rate, the analytic model (Equation 9) calculates the drag force from gas flow properties at infinity. However, the subgrid model implemented in RAMSES evaluates the drag force on the sink based on the same local massweighted quantities used for BHL accretion, and .
In this paper, we thus distinguish between four different drag force calculations:

is the drag force given by Equation 9, using the quantities at infinity.

is the drag force given by Equation 9, using massweighted quantities from the accretion region, i.e. the drag force as estimated by the subgrid model.

is the net gravitational force exerted on the sink by gas on the grid. It is found by summing over all cells in the box, where and are the cell mass and cell position vector relative to the sink respectively. As the setup is symmetric, this is equal to the total gravitational attraction of the wake downstream of the sink.

is the total drag force acting on the sink when the subgrid algorithm is active.
2.5 Plotting conventions
Streamlines annotated on slice plots follow the velocity field in the plane of the slice. They therefore do not represent three dimensional particle trajectories and are merely added to guide the eye as to the dominant flow patterns in the plot. For clarity, a minimum distance is enforced between streamlines. This discrete sampling can mean that some streamlines are interrupted.
3 Accretion
The BHL formula in Equation 3 is based on two limiting cases. The pure Bondi problem, where the sink is at rest relative to the uniform density background, i.e. , and the HoyleLyttleton problem, where . We investigate both cases in detail, in this section and the next, before exploring the parameter space further.
3.1 The Bondi problem
As expected from the analytic work by Bondi & Hoyle (1944); Bondi (1952), Figure 3 shows that the global flow pattern is symmetric and directed radially towards the sink even at very low resolutions, with the sink located at the density peak. Only the immediate sink environment is influenced by the resolution. There is an initial period during which the simulation is dominated by initial condition transients, which last until (see Figure 4). After the simulation has settled into a steady flow pattern, the accretion rate onto the sink, , converges. Note that the actual accretion rate onto the sink, is not necessarily the same as the BHL accretion rate evaluated from the local gas properties , as is limited by the total gas mass present in the accretion region, .
Figure 5 demonstrates that for all simulations presented here, the radial density, pressure, temperature and velocity profiles provide excellent agreement with the analytic predictions by Bondi (1952) at radii larger than the size of the accretion region, i.e. at . Within the accretion region, density and pressure drop sharply as mass is removed from the grid and added to the black hole, while temperatures can increase by up to a factor of 2 above the analytic value. As the analytic profiles diverge as , some deviation within the central region is to be expected, particularly in the velocity field which has to go to zero at the centre of a spherically symmetric flow.
The analytic solution for the BHL accretion rate assumes that gas properties are measured far from the influence of the sink, a situation that is best captured in the simulations presented here when the cell size is significantly larger than the sink’s gravitational scale radius, i.e. . Figure 4 shows that for the unresolved case, here m0.0n0.05a, the accretion rate onto the black hole, as well as the local gas properties as calculated using the cloud particles, reflect the analytic solution closely. The local overdensity before accretion is very shallow, at (see m0.0n0.05a Figure 3, top left), and is suppressed within the cell containing the sink once accretion starts.
The other extreme is the most highly resolved case probed here (m0.0n131a, bottom right panel of Figure 3), where an overdense peak develops within the Bondi radius. When the gas crosses the sonic radius, it transitions from subsonic to supersonic flow as it evolves towards a free fall solution before being accreted. The rising local density increases the BHL rate computed on the fly and the maximum amount of gas permitted is removed from the central cells at each time step. This means that the accretion algorithm effectively transitions from the BHL algorithm to a supply limited accretion scheme, where the accretion rate onto the sink is set by the gas inflow rate into the spherical accretion region of the sink, which has a radius of . We refer to this as Supply Limited Accretion (SLA) throughout this paper.
Note that at this resolution, the accretion rate onto the sink settles well below the analytic BHL rate, to , in agreement with results in Edgar (2004). Simulations with better force resolution probe the density profile on smaller scales, and therefore measure higher densities . As the contraction is adiabatic, higher densities have correspondingly higher sounds speed , as can be seen in Figure 4. For this reason, although the analytic model postulates that the gas should transition to supersonic near the Bondi radius, the black contours in Figure 3 show that the transition occurs at a smaller radius, again in agreement with previous numerical simulations. Finally, the accretion rate has already converged to its value even at the comparatively modest resolution of , in simulation m0.0n32a.
The partially resolved case, m0.0n1.6a, where shows intermediate behaviour, with a shallower central density feature and a smaller evacuated region. This simulation also shows more noticeable grid effects, both in the stream lines and in the central density peak, as spherical symmetry is poorly described by the small number of Cartesian resolution elements when . The resulting steady state accretion rate is lower than the converged value, with as the local density feature feeding the black hole is not replenished efficiently.
Arguably the most worrisome numerical aspect of transitioning to SLA is that the force due to the pressure gradient artificially created by the low density region developing in the immediate vicinity of the sink might dominate the gravitational force on the gas. Figure 6 shows that while such a pressure gradient does reinforce the gravitational pull on the gas at the edge of the accretion region, it is not the dominant force. Morevover, the contribution of this pressure force decreases for simulations with higher resolution as the gravitational acceleration increases faster than the pressure gradient, , within the Bondi radius. As previously mentioned, this effect is also insensitive to our choice of density floor, as long as , as for sufficiently small pressure inside the accretion region, , and therefore only depends on the cell size and the pressure at the edge of the accretion region, , not the actual pressure inside the accretion region .
From the pure Bondi problem presented here we conclude that the accretion algorithm is well behaved at all resolutions. Indeed, at low resolution, when , local gas properties measured in the vicinity of the sink, using cloud particles, produce an accretion rate in excellent agreement with the analytic BHL formula. At high resolution, where , accretion is driven by the local gas supply into the accretion region, and the algorithm transitions to SLA. The accretion rate onto the sink converges to the correct value of in that case, fed by supersonically freefalling gas well within the Bondi radius. For intermediate resolutions of , grid effects lead to poorer spherical symmetry and lower accretion rates, as gas neither reflects the values at infinity nor forms a central gas profile able to efficiently feed the sink, a conclusion also reached by Krumholz et al. (2004). However, independently of resolution, the accretion rate onto the black hole estimated using the subgrid algorithm remains within of the correct value.
3.2 The HoyleLyttleton problem
The other analytic solution was developed by (Hoyle & Lyttleton, 1939) for an accretor moving supersonically through a uniform medium, where the bulk velocity dominates over the local sound speed, such that approaches . In this section we investigate accretion onto the sink in the highly supersonic case where .
3.2.1 The adiabatic case
For the adiabatic case, when , in agreement with expectations from the analytic solutions (Hoyle & Lyttleton, 1939; Ostriker, 1999), we find that a conical wake develops downstream of the sink, as is evident in Figure 7. In the unresolved case, m10n0.08a, the overdensity around the sink is small, and the streamlines and properties of the gas are only mildly perturbed by the presence of the sink. Therefore, , and the sink accretes according to the analytic solution (top panel in Figure 8). The gravitational wake is especially prominent because the simulations presented here are isolated. If the sink was embedded in a nonuniform medium, such as is typically found for black holes in galaxy simulations, we expect that local inhomogeneities would quickly wash out the gravitational focusing effect of the black hole.
With increasing resolution, such as m10n13a in Figure 7, the flow patterns resembles the analytic solution by Hoyle & Lyttleton (1939), with bent streamlines, stagnation point, and accretion column clearly visible in the density slice. The bowshock with the characteristic increase in density towards the edge of the shock, predicted by Ostriker (1999), also becomes apparent. At these intermediate resolutions, the shock is attached to the accretor, and the solution is stable (see Figure 8). As in the pure Bondi case, , and the accretion algorithm transitions to SLA, with the underdense accretion region visible in the density slices (Figure 7).
With even more resolution, and an even smaller accretor, such as in m10n26a, the bowshock detaches from the accretor, and the solution becomes unstable, in agreement with results by Ruffert (1995b). In this regime, the wake alternates between episodes when instabilities die down and the flow returns to a more symmetric configuration (top row, Figure 9), and episodes when instabilities dominate flow patterns, which are severely disruptive (bottom row, Figure 9). While the accretion rate begins to fluctuate due to local instabilities, the timeaveraged value remains similar to the case before the bowshock detached (compare m10n13a and m10n26a in Figure 8), as the accretion column feeding the black hole, characterised by a low vorticity region along the axis of symmetry of the wake continues to reform between unstable episodes (see top right hand panel in Figure 9). When not disrupted by instabilities, the gravitational attraction of the sink creates a peaked density profile around the accretion region, which, like in the Bondi case previously discussed, leads to a resolution dependence of the sound speed at the edge of the accretion region.
While the global wake remains stable, eddies develop behind the shock, suggesting that the instability is caused by the physical acousticadvection instability reported in Foglizzo et al. (2005) for , and “sufficiently small accretors”. According to these authors, the instability is caused by entropy perturbations that advect from the shock to the sonic surface around the accretor, where they excite acoustic waves due to the inhomogeneity of the flow, which in turn propagate outwards back towards the shock, where they excite new entropy perturbations. We would therefore expect the instabilities to occur in the subsonic region, where the gas has been decelerated by the shock and not yet sufficiently reaccelerated by the gravitational potential of the sink.
In the simulations presented here, entropy perturbations (yellow contours) form in the subsonic region between the accretor and the bowshock (see m10n26a and m10n209a in Figure 10) for unstable simulations, as expected from the theory by Foglizzo et al. (2005). For stable simulations, entropy perturbations are concentrated downstream of the accretor instead and are much smaller in magnitude (m10n13a in Figure 10). The origin of instabilities is explored further in Section 5.
The accretion algorithm creates a low density region around the sink, surrounded by a high density shell replenished by inflowing material. While most of the mass that enters the accretion region is removed, the sink particle algorithm does not implement a strict inflow criteria, with any surplus gas free to leave the accretion region during the next timestep. Accretion onto the sink varies on short timescales, as the turbulent flow feeds the black hole intermittently (see m10n26a in Figure 8), with an average value of , in agreement with simulations discussed in Edgar (2004). This is also the value recovered for the steady state solution of extremely small accretors () in the 2D axisymmetric simulations of El Mellah & Casse (2015).
The picture changes drastically for the highest resolutions, and smallest accretors, (m10n209a in Figure 7), as the instabilities become more severe. The bow shock opening angle increases, supported by a strong rotational flow around the sink, with vorticity near the sink increasing by up to an order of magnitude (see Figure 10). The accretion column, which feeds the sink at lower resolutions, is permanently disrupted and the solution is entirely dominated by instabilities, as evident by the oscillations in Figure 8. Due to the increasing vorticity behind the bow shock, the gas within builds up coherent angular momentum (top panel, Figure 11). This angular momentum provides nonnegligible rotational support against gravitational collapse for m10n209a (the spin parameter, , where and are the total mass and angular momentum contained within a sphere of radius , and is the circular velocity at a given radius (Bullock et al., 2001), bottom panel, Figure 11), preventing gas from accreting efficiently. Periods of lower angular momentum around the black hole, seen for resolutions , disappear. As a result, the time averaged accretion rate drops to , with order of magnitude fluctuations around this value, so that extrapolating results from lower resolution runs provides a poor estimate of the accretion rate onto the sink.
These results are in agreement with work by MacLeod & RamirezRuiz (2015), who investigate HoyleLyttleton type accretion in the presence of a density gradient at infinity, and find that the resulting circularisation of gas behind the (warped) shock reduces the accretion rate onto the sink by over an order of magnitude. In the simulations presented here, the circularisation is driven by the advectiveacoustic instability, not a global density gradient. However, both cases shows that the presence of significant angular momentum in the gas behind the shock reduces the accretion onto the sink by more than order of magnitude.
Overall, we conclude that, in contrast with the Bondi case, the accretion onto the sink for an adiabatic HoyleLyttleton flow seems to converge to a value well below the analytic BHL solution, as it becomes entirely dominated by instabilities for very small accretors, i.e. (equivalent to ).
3.2.2 The quasiisothermal case
Foglizzo et al. (2005) argue that the advectiveacoustic instability should disappear for , whereas Ruffert (1995a) report unstable flow for at any . To further investigate the possible presence of instabilities, we ran a suite of quasiisothermal simulations with . We find a gas flow pattern broadly similar to the adiabatic case presented above, albeit with very different consequences for the accretion onto the sink. As can be seen in Figure 13, at very low resolution, for m10n0.08i, only a small overdensity forms downstream of the sink, which accretes according to the analytic BHL formula (see the top panel in Figure 12). For higher resolutions, , the local flow pattern again shows curved streamlines, stagnation point, shock and SLA. However, the shock opening angle is significantly smaller. Figure 12 also shows that the accretion onto the sink is unstable for , a significantly lower resolution threshold than for the runs. Like in the adiabatic case, simulations with an intermediate resolution (e.g. m10n3.3i) have a reduced accretion rate as the resolution does not yet allow the characteristic flow pattern to fully emerge. As is discussed further in Section 5, at a resolution of , the stagnation point lies within the accretion region of the accretor, slowing down gas flows onto the sink. By definition, the sound speed (bottom panel of Figure 12) remains constant in the quasiisothermal case, so as resolution increases, the relative velocity of the gas with respect to the sink increases much more than in the adiabatic case (bottom panel of Figures 8 and 12).
As a result, and contrary to the adiabatic case, the shock remains attached to the accretion region regardless of resolution and the flow stays supersonic everywhere except in a small narrow region around the stagnation point (see Figure 13). Despite this, we find that strong instabilities develop in the wake, particularly for small accretors, which also lead to accretion rate variations on short timescales. Instead of originating at the bowshock, these instabilities occur downstream of the sink and affect the wake globally, disrupting the accretion column onto the sink. No eddies are visible in the stream lines (see Figure 13). Rather, the narrow wake clumps and is distorted in the direction perpendicular to the axis of symmetry of the wake, in agreement with the seminal work of Ruffert (1995b). Taken at face value, the persistent instabilities for seem to contradict predictions by Foglizzo et al. (2005), who argue that the acousticadvective instability should disappear for . However these authors do point out that a different type of instability could occur at low , based on vorticity (rather than entropy) perturbations between shock front and accretor. We note that the instabilities seen here are also reminiscent of the ones discussed in Cowie (1977) and Soker (1990), who study the fate of small overdensities when the wake is modelled as an accretion line, but postpone a further discussion of the origin of these instabilities to Section 5. Beyond the exact nature of the instabilities, the main difference between the adiabatic and the quasiisothermal simulations is that the strong drop in average accretion rate for the smallest accretors, which have , is absent for quasiisothermal simulations. At the highest resolution probed here, the accretion rate converges to : while instabilities exist, they do not efficiently prevent accretion onto the sink.
In summary, we conclude that for black hole moving supersonically, with , accretion proceeds via the BHL algorithm at low resolution, , and transition to SLA at higher resolution. The wake is unstable for both values of , with the instabilities in the adiabatic case originating in the subsonic region between the sink and the detached bowshock, and near the stagnation point in the quasiisothermal case. With increasing resolution, and decreasing size of the accretor, the flow becomes unstable on progressively shorter timescales. In the adiabatic case, instabilities dominate for , which leads to an order of magnitude reduction in the timeaveraged accretion rate onto the sink, whilst in the quasiisothermal case this averaged accretion rate converges to a value close to that of the BHL case.
3.3 Exploring the full parameter space in
We now explore the evolution of flow patterns and accretion rates in the adiabatic case for a wider range of Mach numbers and different resolutions. Note that in the case of intermediate Mach numbers (), we do not expect the analytic BHL formula (Equation 3) to yield as accurate an estimate of the accretion onto the sink as in the low and high Mach number cases previously studied as it is merely an educated interpolation between these two extremes cases.
Having said that, Figure 15 presents the average accretion rate, in units of the analytic BHL accretion rate , for a variety of Mach numbers, as a function of resolution. As could already be seen for the Bondi and the HoyleLyttleton problem (Section 3.1 and 3.2 respectively), the unresolved case with closely follows the analytic formula. As resolution increases, the behaviour becomes more complicated. In the intermediate regime, where , the simulations diverge from the analytic formula in a way that depends nonmonotonically on the Mach number. Sub and supersonic simulations ( and respectively) show average accretion rates systematically lower than the BHL formula, by up to a factor . Transsonic simulations, on the other hand, with feature accretion rates which are larger by up to a factor of . The BHL formula, used here to normalise results, is most uncertain for the transsonic regime, where both the bulk velocity and the sound speed have a significant influence on the flow. Our high resolution results () support this conclusion as the accretion rates indeed converges to higher values in the transsonic regime. We caution that at intermediate resolutions, the pressure force into the low resolution region around the sink can dominate over the local gravitational force, possibly funnelling extra gas into the accretion region and thus leading to an overestimate of the accretion rate. However, this effect is very localised as it only occurs at the edge of the accretion region, and is alleviated by the kernel function which creates a gradual transition of density within the region covered by the cloud particles. However at high resolution (), the gravitational force comfortably dominates over the pressure force for all cases studied here (see Figure 6 for a measure in the Bondi case). Moreover, all resolved simulations with Mach numbers show steady state solutions (Figure 14).
From on, eddies begin to form behind the shock, and when the accretor becomes small enough, instabilities becomes stronger and begin to influence accretion onto the sink more significantly. The density slices for both and show strong instabilities that disrupt the flow patterns and decrease the timeaveraged accretion rate onto the sink by up to an order of magnitude below the analytic value. We caution that while a lot of care has been taken to minimise the impact of initial conditions (see Appendix A), the seeding of the instabilities could still be due to the way the simulations are initialised. This is likely to affect the exact resolution and/or Mach number at which the instability dominated regime appears, but unlikely to make it vanish altogether.
In summary, we conclude that the sink particle algorithm, using a locally evaluated BHL accretion rate as described in Section 3, is a versatile subgrid model that smoothly adapts to a variety of resolutions. For highly resolved simulations (), the kernel function used to remove mass ensures that the maximum accreted mass per timestep, dominated by the dense cells at the edge of the accretion region, always exceeds the local gas supply and the subgrid model automatically transitions to SLA. Intermediate and low resolutions () lead to mixed results and appear to be a difficult regime when approximating the accretion rate onto the black hole from local gas properties, but still manage to capture the accretion rate within a factor , at least in the moderate Mach number regime (). It is only at higher Mach numbers that they deviate from resolved time averages by more than an order of magnitude, an effect that can significantly impact the final mass of the sink. This is potentially important when simulating the cosmological growth of supermassive black holes, where an early accretion boost is crucial to enable the black holes to reach observed masses within the limited timeframe available (See Volonteri, 2010, for a review).
4 Drag force
The gravitational force of the wake formed downstream of sink particles in the presence of a significant bulk flow exerts dynamical friction on the sink particle, opposite to the direction of motion, causing it to reduce its relative velocity with respect to the gas over time. However, the total drag force of the wake depends sensitively on the mass contained in the wake, particularly close to the sink.
Returning to the adiabatic simulations discussed in detail in Section 3.2.1, one can see in the density slices in Figure 7 that the wake develops even for low resolution simulations. Figure 16 shows the mass distribution and drag force profile of each wake at , plotted against the distance to the sink measured along the axis of symmetry of the wake. The mass distribution close to the sink, where , shows some variation with increasing resolution, due to the locally unstable flow. However, the global structure of the wake at larger radii converges quickly, for , and remains stable. The highest contributions to the drag force are found in the immediate vicinity of the accretor, but as symmetric contributions upstream and downstream of the sink cancel out, the larger scale structure of the wake (within in this case) contributes the bulk of the net drag force onto the sink. As a result, the overall drag is adequately captured even with moderate resolution: integrating over for simulations m1013a, m10n26a, m10n52a, m10n209a yields values for of 0.057, 0.063, 0.071, 0.059 respectively in the dimensionless units used in this work, so that even in the most violently unstable case (m10n209a), the drag force fluctuates by less than 20%.
Figure 17 shows the spatial contribution of density slices along the wake for adiabatic simulations at for a range of Mach numbers, the density slices for which can be found in Figure 14. As expected, the wakes contain a significant amount of mass on relatively large scales (up to ), especially for transsonic () configurations. However, the inverse square dependence on the distance to the sink means that for all Mach numbers investigated here, most of the gravitational drag force comes from a region within . Note that the intensity of the force also depends on the opening angle of the wake, with a similar mass profile exerting a stronger pull in the direction of motion if confined to a narrower cone. Finally, we also measure a nonnegligible contribution of material in front of the sink particle, pooling behind the detached shock, that exerts a gravitational force in the opposite direction and reduces the overall drag, especially in the sub and transsonic regimes. While this feature was also observed in Chapon et al. (2013), it appears more prominently in the simulations presented here, and is completely absent in the analytic solutions for supersonic black holes by Ostriker (1999), who state that the sink particle only generates a density wake within the rear Mach cone.
Comparing the total net gravitational drag for the set of resolved simulations in Figure 17 to analytic estimates in Figure 18, and considering that the size of the accretor sets the smallest scale , the Coulomb logarithm evaluates to , larger than the value of reported by Chapon et al. (2013) in the transsonic regime. However it is clear from the bottom panel of Figure 17 that , the characteristic size of the medium which the accreting object traverses, is quite a sensitive function of Mach number, as it drops from in the transsonic regime to a value of about at Mach , corresponding to . Fitting over the whole range of Mach numbers covered by our simulations, we find a time averaged best fit Coulomb logarithm of , with (5.0) providing an adequate lower (upper) bound to individual time measurements. We suspect that this discrepancy with the Chapon et al. (2013) results is partly due to the different value of used in these authors’ simulations, but more likely caused by the absence of accretion onto their black holes. This is somewhat corroborated by the fact that we also find an excess in gravitational drag for subsonic sinks, which show a more prominently asymmetric density profile than in the nonaccreting analytic solutions.
Contrary to the accretion rates in Figure 15, the drag force due to the wake forming behind the sink shows rapid convergence at a surprisingly low resolution of , as shown in Figure 19. It is only vastly underestimated at resolutions as low as , i.e. when the scale length corresponding to the gravitational influence of the black hole, , is small in comparison to the minimal cell size, so that gravitational focusing is inefficient. The drag force due to the wake is only moderately influenced (up to ) by the instabilities developing behind the bow shock, as it is dominated by larger scale contributions.
On top of the direct computation of the gravitational drag exerted on the sink by the overdensity in its wake, RAMSES includes a subgrid algorithm that calculates the drag force based on the same local mass weighted average quantities used to estimate the accretion rate, , with . It is designed to compensate for the lack of drag force in low resolution simulations, clearly visible in Figure 19 for . As the drag force on the sink is a vector quantity, it inherits its direction from the local relative velocity, .
For the supersonic case, , the top panel of Figure 20 shows that the magnitude of the cell based force, converges for and . The instabilities in the flow for cause small variations of about in magnitude (filled symbols). To investigate the impact of resolution both on the magnitude and direction of the drag force, we define a fiducial drag , based on the flow velocity at infinity and the converged drag force intensity . Figure 20 shows that the drag force due to the wake is always parallel to the fiducial force (filled markers, middle panel). As expected from the direction of the wake, the force acts in the opposite direction to the global flow velocity, and slows the sink down. For low resolution, , the drag force due to the wake can be underestimated by a large factor, but remains steady. For high resolutions, , the magnitude of the force rapidly converges to within of the fiducial value. At the highest resolutions, , short term variations in the magnitude of the force are visible (filled markers, top panel) but no significant deviation from the axis of symmetry of the wake (middle panel).
The subgrid based drag force, (see top two panels of Figure 20), by contrast, shows a very erratic behaviour for . Its magnitude, (solid lines), varies considerably on the shortest timescale probed here, the finest timestep of the simulations, and can be both significantly larger or significantly smaller than the value. Moreover, higher resolutions show larger fluctuations. Not only does the magnitude of the force fluctuate rapidly, but it is also directed in the opposite direction to the fiducial force most of the time (actually at all times for ). This can be easily understood, because the bulk of the mass enters the accretor through the accretion column, which has a flow direction directly opposite to the global flow (see Figure 7 for some examples). As a result of mass weighting, the velocity of the accretion column thus dominates the local flow velocity, , and as , the subgrid drag force flips direction as soon as the accretion column forms. The full extent of the problem becomes apparent when calculating the total drag force in the presence of the subgrid algorithm, (bottom panel of Figure 20). For the resolved cases, , frequently is the dominant term, causing a net force that accelerates the sink relative to the global gas flow. This is clearly unphysical, and entirely caused by the fact that the local mass weighted relative velocity, does not reflect either the direction and/or magnitude of the value at infinity, as soon as the accretion column forms and the bow shock detaches from the accretor. As expected, the total drag force for the unresolved cases is unaffected, as the local velocity reflects the value at infinity because the accretion column has not (fully) developed. The subgrid based drag force therefore significantly and accurately contributes to the overall drag force in that case. This contribution naturally drops as the accretion column builds up, leading to a decrease in relative velocity and an increase in sound speed as the flow in the vicinity of the accretor becomes subsonic (see bottom panel of Figure 8). From Figure 20 (top and middle panel) the flip in drag force direction occurs around at , which does not create a significant problem at this resolution as the resolved drag force term dominates. As a result, we find (bottom panel of Figure 20) for the standard drag force subgrid model used in RAMSES with , and near perfect agreement at with our best fit value of (not shown).
Based on these results, and taking a conservative approach, subgrid algorithms for the drag force onto the sink should be avoided as soon as the characteristic scale radius, , becomes larger than the size of the accretor, . Note that itself depends on the relative velocity of the black hole and the ISM and can therefore be difficult to determine in more complex simulations than the idealised BondiHoyleLyttleton flows investigated here, where the values at infinity are known. The criteria at which the subgrid drag force on black holes becomes unphysical might therefore have to be revised for galaxy evolution simulations, the analysis of which we leave for future work.
5 Are the instabilities at high Mach number physical?
To the best of our knowledge, no complete analytic analysis yet exists to explain the instabilities regularly found for accretion wakes in the HoyleLyttleton problem, so the question remains whether the observed instabilities are physical or numerical in nature. Our work on the danger of using small accretors with uniform initial conditions (see Appendix A) shows that a numerical origin is extremely difficult to rule out. However, a vast amount of efforts to stabilise simulations have failed, as instabilities have been reported using a wide variety of codes, coordinate systems and models for the accretor (See Foglizzo et al., 2005, for a review). In agreement with the results presented here, a majority of authors report a clear link of the appearance of instabilities with the size of the accretor.
Claims of the existence of steady state solutions for small accretors and high Mach numbers, such as the work by El Mellah & Casse (2015) ( and ) and Pogorelov et al. (2000) ( and ) have also been made. However, all such results are based on 2D axisymmetric simulations, and it is well known that for HoyleLyttleton simulations, certain instabilities only appear for particular configurations. For example the flipflop instability, frequently observed in 2D, is entirely absent in 3D simulations (Blondin & Pope, 2009).
To investigate how a shrinking accretor impacts the stability of the flow, we take an established, stable simulation with an intermediate accretor size (m3n32a), and shrink the accretor by adding refinement levels to reproduce the conditions in m3n72a, where instabilities unsettle the flow. Both simulations are supersonic, with , and m3n32a has reached a steady state before the accretor is shrunk.
The density slices in Figure 21 confirm observations from the case presented in Section 3.2. Shrinking the accretor creates a subsonic region, forming upstream between the supersonic region surrounding the accretor and the shock front, which is absent for the steady state solution at lower resolution, as can be seen in Figure 21. As soon as the subsonic region develops, entropy perturbations form close to the accretor upstream of the sink (Figure 22). They oscillate back and forth between the two sonic surfaces, disturbing the upstream bowshock in the process. As can be seen both in the profiles in Figure 22 and in the contours in Figure 23, the entropy perturbations remain confined between the two sonic surfaces, and therefore are only able to form if the accretor is sufficiently small to allow this subsonic region to open. During the next 10 dynamical times, the simulation does not resettle into a steady state.
As previously mentioned, Foglizzo (2009) suggest that the advectiveacoustic instability should unsettle the flow in exactly this manner. It is characterised by small entropy perturbations forming behind the shock, which are advected towards the accretor, where the local rise in density causes them to reflect back to the shock front, exciting further perturbations. Contrary to Foglizzo (2009), who predict that entropy perturbations form at the bowshock and propagate towards the accretor, the simulations presented here show that the entropy perturbations first appear close to the accretor and then propagate forwards in the direction of the bow shock instead. As was shown in section 3.3, all of our simulations which have , a sufficiently strong shock (which occurs for ), and enough resolution for the shock to detach, show these types of instabilities, providing strong support for the physical origin advocated in the analysis of these authors.
While the entropy perturbations upstream of the sink appear as soon as the subsonic region opens up, the flow downstream of the sinks takes longer to become unsettled. This can be seen in Figure 24, which shows the minimum vorticity within a hemisphere of radius , centred on the sink and oriented downstream, plotted against the instantaneous accretion rate of the sink particle. On timescales of order 5 flow crossing times of , the flow at the bottom end of the accretion column evolves from a lowvorticity state, with , in which the black hole is accreting efficiently at to a high vorticity state with , in which accretion rates vary strongly and can drop as low as (Figure 24). While the vorticity perturbations behind the bowshock are difficult to quantify due to the nonuniform, resolution dependent vorticity pattern of the characteristic HoyleLyttleton flow, the vorticity slices in Figure 23 show that both vorticity troughs and peaks develop in the regions affected by entropy perturbations (yellows contours) as the flow forms drawn out eddies (see flow lines in Figure 21). Over time, these eddies stretch backwards along the wake, eventually upsetting the symmetry of the flow downstream of the accretor, including the accretion column feeding the sink.
For the isothermal case, the shock remains attached and the flow is supersonic everywhere, but the wake again becomes globally unstable for sufficiently small accretors. This is particularly important for simulations of SMBHs, as their galactic environment is strongly affected by radiative cooling. While the instabilities for the quasiisothermal case appear to have a different — and arguably less well analytically understood — origin than the adiabatic ones, they lead us to conclude that Hoyle–Lyttleton accretion onto galactic black holes is very unlikely to reach a steady state, whether radiative cooling plays an important role or not.
Indeed, repeating our experiment to explore the origin of the instability in the quasiisothermal case, by shrinking the accretor for a simulation from (m3n3.5i), which leads to a steady state solution, to (m3n7i), which is unstable (see Figure 25) reveals that the instability originates downstream of the sink in the region of the stagnation point of the flow. Instead of an increase in entropy in front of the accretor, like in the adiabatic case, the quasiisothermal case shows a low entropy region within the entire bowshock (Figure 26). This low entropy region is present for both stable and unstable simulations, with no signs of entropy perturbations, and we therefore conclude that this instability is not linked to the advectiveacoustic cycle found in the adiabatic case.
Instead, this could be one of the two instabilities of the accretion column discussed in Foglizzo et al. (2005): either the longitudinal instability first described in Cowie (1977), where the authors consider the impact of a small density perturbation in the accretion column, or the transverse instability of the accretion column from Soker (1990). Investigating the longitudinal case first, Figure 27 shows the density contrast along the wake. As predicted by Cowie (1977), density perturbations form near the stagnation point but we do not see the amplification of waves travelling towards the accretor discussed by Foglizzo et al. (2005). Instead, density perturbations propagate downstream from the accretion point, and have all but faded by , with no evidence of a selfsustaining feedback loop. We therefore conclude that the longitudinal instability is not responsible for unsettling the wake.
Figure 28 shows that the stagnation point also plays an important role in the transverse displacement of the wake. The threedimensional angle , measured between the axis of symmetry of the problem and centre of mass of a given slice along the wake, is largest in the region bounded by the stagnation point and the accretor. This confirms the visual conclusions from Figure 13, where the largest perpendicular displacement of the wake occurs at or near the stagnation point for all three simulations with . We see no evidence for the predicted increase of amplifications with increasing distance from the accretor (Soker, 1990), and in fact see little transverse displacement of the wake beyond the stagnation point at all. We therefore conclude that while a local transverse instability can be observed in the threedimensional simulations presented here, its effect is constrained to the supersonic region between the stagnation point and the accretor. If this region is unresolved, such as in cases of very low resolution, the accretion column, and by extension the wake, remain stable.
The accretion column also remains stable in the adiabatic case previously discussed, as both the longitudinal and the transverse instability rely on the assumption that the radius of the shock cone is small compared to other distances, i.e. that the wake is very narrow. While theoretically both instabilities should appear for sufficiently high Mach numbers in the adiabatic case, the advectiveacoustic instability acts to increase the shock opening angle for moderate to high Mach numbers, suppressing the instabilities of the accretion column.
Contrary to the adiabatic case, accretion onto the black hole remains effective even in the presence of turbulence. As can be seen in Figure 29, gas properties in the wake remain in the low vorticity  high accretion regime, where , throughout the course of the simulation despite the fact that the isothermal case was run for an extra 10 dynamical times. The vorticity slices in Figure 26 confirm that a low vorticity region continues to exist within the narrow bow shock, allowing gas to continue feeding the black hole.
A variety of other instabilities have been suggested to potentially arise from HoyleLyttleton ’like’ accretion flows, such as individual vortex rings in Kim & Kim (2009) for a nonaccreting massive perturber embedded in an adiabatic 2D supersonic flow, or shock cone vibrations as reported in LoraClavijo & Guzman (2013) when investigating 2D adiabatic BHL accretion using general relativity instead of Newtonian gravity. However, neither of these phenomena are observed here. Finally, we stress that the main numerical limitation of the simulations presented here is that the size of the accretor is intrinsically tied to the resolution, with . To truly investigate the nature of the instabilities would require exploring the impact of resolution and accretor size independently, which is beyond the scope of this work. We expect that this shortcoming will only affect the details of the instability triggering and their exact intensity.
6 Conclusions
As numerical capabilities increase, we are able to cover a larger range of length scales in simulations. In this paper, we explored the impact of resolution on the accretion and drag force for a sink particle in grid simulations, using subgrid algorithms based on the BondiHoyleLyttleton (BHL) interpolation formula for accretion, and the linear analytic drag force estimate of Ostriker (1999). In all simulations presented here, the size of the accretion region shrunk with increasing resolution, and both analytic formulas were evaluated using massaveraged, kernelweighted quantities in the immediate vicinity of the accretor. We adopted an idealised BHL setup, where originally homogeneously distributed gas with a uniform Mach number was allowed to settle in the gravitational potential of the accreting particle, with the relevant scale radius, , resolved by a fixed number of resolution elements, throughout the simulation. We then investigated the impact of a wide range of resolutions, , Mach numbers , and two values of the adiabatic index , on the accretion rate and drag force subgrid algorithms onto our accretors, modelled as sink particles with radius of .
We found that, as expected, for very low resolutions, i.e. , where the local mass weighted quantities reflect the ones at infinity because the gravitational influence of the sink is small, the accretion subgrid algorithm closely followed the analytic values. Accretion rates converged for and Mach numbers , although not to the analytic value set by the accretion algorithm as accretion onto the black hole transitioned from BHL to supply limited accretion. However, in this situation, the accretion rates measured onto the accretor all converged to values within a factor of the BHL analytic interpolation formula, with the largest discrepancies found for transsonic flows.
On the other hand, at higher Mach numbers, , and higher resolution, instabilities appeared and started to dominate the flow to such an extent that the accretion rate dropped to in the adiabatic runs with . In the quasi isothermal case (), accretion rates remained within a factor of the BHL rate over the whole range of Mach numbers probed here, despite the presence of instabilities in the accretion column at high Mach numbers ().
These findings are in agreement with previous work by Ruffert (1995b); Foglizzo et al. (2005) as to the origin of these instabilities, at least in the adiabatic case when . They are caused by an acousticadvective feedback loop which develops in the subsonic region between the shock front and the accretor, as described analytically in Foglizzo et al. (2005); Foglizzo (2009). In the quasiisothermal case, where , we instead observe a transverse instability of the accretion column in the region between the stagnation point and the accretor, which locally displaces this column perpendicularly to the axis of symmetry of the problem, reminiscent of the transverse instability discussed in Soker (1990).
The drag force due to the wake, calculated from the cell density in the box, converged quickly with resolution, again in agreement with results by Ruffert (1995a, b). However, as the bow shock and accretion column began to be resolved by the simulation, the local mass weighted quantities used in the subgrid algorithm to evaluate the analytic drag force according to the formula of Ostriker (1999) started to poorly reflect the values of these quantities “at infinity”. Indeed, not only did the magnitude of the drag force estimated in this way fluctuate on extremely short time scales, but its direction actually flipped, so that even after adding the contribution from the resolved wake, the net force on the accretor remained an acceleration with respect to the gas flow, rather than a drag. This clearly unphysical behaviour solely occurs because the local relative velocity becomes dominated by gas flowing down the accretion column. Therefore, subgrid algorithms for the drag force should be avoided as soon as the characteristic radius becomes larger than the size of the accretor.
The main limitation of the conclusions we reach about the subgrid algorithms presented in this work, is that quantities “at infinity” are poorly defined in more realistic galaxy evolution or cosmological simulations. For example, we have measured the scale radius for quantities “at infinity”, taken to be the boundary of the box. However, embedded in a galaxy, the flow far from the accretor is not uniform or isothermal, and the gravitational field due to the accreting black hole not necessarily dominant. We also caution that it is important to resolve , which can be much smaller than the Bondi radius for highly supersonic flows which are quite frequent in galaxy simulations.
Moreover, while we measure fairly similar trends in both the adiabatic and quasiisothermal cases probed here, our simulations do not include radiative cooling explicitly, which could significantly influence the sound speed in the vicinity of the accretor, and through it both the accretion rate and drag force. In addition, while we calculate the drag force, we do not allow the sink to move under its influence. This is particularly important for the instability dominated runs, where we would not expect the accretor to remain located on the axis of symmetry of the bow shock at all times. In less idealised simulations than those presented in this work, we also expect the sink to encounter nonhomogeneous gas moving nonuniformly, other compact objects as well as gradients in the local gravitational potential.
However, it is clear that the subgrid models we have described in this work go beyond a basic implementation of an analytic BHL prescription. They are more versatile and reliable, capable of handling a wide range of resolutions and complex flow configurations. In particular, the smooth transition from an algorithm based on the analytic BHL accretion formula to SLA when local density and velocity features are resolved ensures that the accretion rate onto the black hole makes the best use of the local information available in the simulation, and should quite naturally converge to the correct solution once the size of the accretion region approaches the physical size of the black hole. By contrast, the drag force subgrid algorithm becomes unphysical as resolution increases and extra caution regarding its implementation in less idealised simulations is needed, a question which we plan to address in a followup work.
Acknowledgements
The authors thank Will Potter for helpful suggestions, and the reviewer for a constructive referee report that significantly improved the depth and quality of this paper. The research of RSB is supported by Science and Technologies Facilities Council (STFC) and by the Centre National de la Recherche Scientifique (CNRS) on grant ANR16CE310011, and the research of AS and JD at Oxford is supported by the Oxford Martin School and Adrian Beecroft. This work is part of the HorizonUK project, which used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk ). This equipment is funded by BIS National EInfrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National EInfrastructure. All visualisations were produced using the ytproject (Turk et al., 2011).
References
 AnglésAlcázar et al. (2015) AnglésAlcázar D., Özel F., Davé R., Katz N., Kollmeier J. A., Oppenheimer B. D., 2015, ApJ, 800, 127
 Bate et al. (1995) Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362
 Biernacki et al. (2017) Biernacki P., Teyssier R., Bleuler A., 2017, MNRAS, 469, 295
 Blondin & Pope (2009) Blondin J. M., Pope T. C., 2009, ApJ, 700, 95
 Bondi (1952) Bondi H., 1952, MNRAS, 112, 195
 Bondi & Hoyle (1944) Bondi H., Hoyle F., 1944, MNRAS, 104, 273
 Booth & Schaye (2009) Booth C. M., Schaye J., 2009, MNRAS, 398, 53
 Bullock et al. (2001) Bullock J. S., Dekel A., Kolatt T. S., Kravtsov A. V., Klypin A. A., Porciani C., Primack J. R., 2001, ApJ, 555, 240
 Chapon et al. (2013) Chapon D., Mayer L., Teyssier R., 2013, MNRAS, 429, 3114
 Cowie (1977) Cowie L. L., 1977, MNRAS, 180, 491
 DeBuhr et al. (2010) DeBuhr J., Quataert E., Ma C.P., 2010, MNRAS, 412, 1341
 Dubois et al. (2010) Dubois Y., Devriendt J., Slyz A., Teyssier R., 2010, MNRAS, 409, 985
 Dubois et al. (2013) Dubois Y., Pichon C., Devriendt J., Silk J., Haehnelt M., Kimm T., Slyz A., 2013, MNRAS, 428, 2885
 Dubois et al. (2014) Dubois Y., Volonteri M., Silk J., 2014, MNRAS, 440, 1590
 Edgar (2004) Edgar R., 2004, New Astron. Rev., 48, 843
 El Mellah & Casse (2015) El Mellah I., Casse F., 2015, MNRAS, 454, 2657
 Elling (2009) Elling V., 2009, Acta Math. Sci., 29, 1647
 Foglizzo (2009) Foglizzo T., 2009, ApJ, 694, 820
 Foglizzo et al. (2005) Foglizzo T., Galletti P., Ruffert M., 2005, A&A, 435, 397
 Hobbs et al. (2011) Hobbs A., Nayakshin S., Power C., King A., 2011, MNRAS, 413, 2633
 Hopkins & Quataert (2011) Hopkins P. F., Quataert E., 2011, MNRAS, 415, 1027
 Hoyle & Lyttleton (1939) Hoyle F., Lyttleton R. a., 1939, Math. Proc. Cambridge Philos. Soc., 35, 405
 Kim & Kim (2009) Kim H., Kim W.T., 2009, ApJ, 703, 1278
 Korol et al. (2016) Korol V., Ciotti L., Pellegrini S., 2016, MNRAS, 460, 1188
 Krumholz et al. (2004) Krumholz M. R., McKee C. F., Klein R. I., 2004, ApJ, 611, 399
 LoraClavijo & Guzman (2013) LoraClavijo F. D., Guzman F. S., 2013, MNRAS, 429, 3144
 MacLeod & RamirezRuiz (2015) MacLeod M., RamirezRuiz E., 2015, ApJ, 803, 41
 Negri & Volonteri (2017) Negri A., Volonteri M., 2017, MNRAS, 467, 3475
 Ostriker (1999) Ostriker E. C., 1999, ApJ, 513, 252
 Peery & Imlay (1988) Peery K., Imlay S., 1988, in 24th Jt. Propuls. Conf.. Reston, Virigina
 Pogorelov et al. (2000) Pogorelov N. V., Ohsugi Y., Matsuda T., 2000, MNRAS, 313, 198
 Power et al. (2011) Power C., Nayakshin S., King A., 2011, MNRAS, 412, 269
 Quirk (1994) Quirk J. J., 1994, Int. J. Numer. Methods Fluids, 18, 555
 RosasGuevara et al. (2015) RosasGuevara Y. M., et al., 2015, MNRAS, 454, 1038
 Ruffert (1995a) Ruffert M., 1995a, A&A, 113, 133
 Ruffert (1995b) Ruffert M., 1995b, A&A, 331, 817
 Sijacki et al. (2007) Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877
 Soker (1990) Soker N., 1990, Astron. J., 99, 1869
 Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
 Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
 Teyssier et al. (2011) Teyssier R., Moore B., Martizzi D., Dubois Y., Mayer L., 2011, MNRAS, 414, 195
 Turk et al. (2011) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, ApJS, 192, 9
 Volonteri (2010) Volonteri M., 2010, Annu. Rev. Astron. Astrophys., 18, 279
 Volonteri et al. (2016) Volonteri M., Dubois Y., Pichon C., Devriendt J., 2016, MNRAS, 460, 2979
Appendix A Uniform initial conditions and small accretors
The simulations presented in this work use uniform initial conditions to simplify and homogeneise the setup of different flow configurations. After an initial period, during which the flow settles, the simulations are assumed to have erased the memory of their initial conditions and evolved to their natural quasisteady state solution. As can be seen in the top row of Figure 30, small accretors can cause shocks during the initial settling phase which generate density perturbations that remain present even after many dynamical time scales. In this Appendix, we present the method used here to minimise this spurious effect.
The first precaution we take is to introduce the analytic gravity field gradually, with while , where . For large accretors (), the simulations show some instability for , visible in the profile plots of Figure 31, particularly when the bowshock detaches, but then settle into stable and converged solutions. However, for sufficiently small accretors, i.e. when , simulations of any Mach number show a shock forming in front of the sink and travelling upstream, driven and supported by eddies (see top row of Figure 30 for an example). This ’circularisation’ of the gas flow breaks the symmetry of the problem, which develops perturbations that fail to settle during more than 10 dynamical times. We only observed this phenomenon for sufficiently small accretors, where the density profile around the sink is more strongly peaked, causing the reverse shock.
In order to minimise the effect of the initial conditions for small accretors, all simulations in this paper with were run using a reduced maximum level of refinement while , which is equivalent to an effective accretor size . This avoids shocking the gas and allows the flow to settle. Around , the extra levels of refinement are added in a staggered fashion to minimise instabilities caused by changing the grid structure: we first decrease the size of the accretor and then add an extra level of refinement per timestep until the target resolution is reached.
Figure 31 shows that some instabilities do occur during and just after adding refinement levels, as evidenced by the high frequency fluctuations in the gas properties , and caused by the gas “sloshing” at the bottom of the gravitational well. After a few dynamical times, the solution settles into the stable flow pattern observed for larger accretors (bottom row, Figure 30). A test run using (m0.5n65a), which is naturally stable, shows that both the run with the full resolution available from the beginning (n65_t0 in Figure 31), and the run where the full resolution is reached at (n65_t3 in Figure 31) converge to , which is also recovered for the small accretor with staggered refinement, n263_t3. By contrast, the simulation in which the small accretor is fully resolved from the start, n262_t0, displays a lack of convergence, high frequency variations in the accretion rate and a much lower time averaged value of . To ensure that simulations have sufficient time to settle, we measure all time averaged values in the paper for .