Outflows from inflows: the nature of Bondi-like accretion

ABSTRACT

The classic Bondi solution remains a common starting point both for studying black hole growth across cosmic time in cosmological simulations and for smaller scale simulations of active galactic nuclei (AGN) feedback. In nature, however, there will be inhomogeneous distributions of rotational velocity and density along the outer radius (R_o) marking the sphere of influence of a black hole. While there have been many studies of how the Bondi solution changes with a prescribed angular momentum boundary condition, they have all assumed a constant density at R_o. In this Letter, we show that a non-uniform density at R_o causes a meridional flow and due to conservation of angular momentum, the Bondi solution qualitatively changes into an inflow–outflow solution. Using physical arguments, we analytically identify the critical logarithmic density gradient ||$\partial \ln \rho/\partial \theta$|| above which this change of the solution occurs. For realistic R_o, this critical gradient is less than 0.01 and tends to 0 as R_o → ∞. We show using numerical simulations that, unlike for solutions with an imposed rotational velocity, the accretion rate for solutions under an inhomogeneous density boundary condition remains constant at nearly the Bondi rate |$\dot{M}_\mathrm{ B}$|⁠, while the outflow rate can greatly exceed |$\dot{M}_\mathrm{ B}$|⁠.

accretion, accretion discs, black hole physics, hydrodynamics

1 INTRODUCTION

There are multiple modes of accretion thought to be important for black hole (BH) growth. The so-called hot accretion mode refers to disc mediated accretion through geometrically thick discs that are either optically thin (as in advection-dominated accretion flows, ADAFs; Narayan & Yi 1994; see Yuan & Narayan 2014 for a review) or optically thick (as in slim discs; Abramowicz et al. 1988). Likewise, the cold accretion mode refers to standard geometrically thin, optically thick discs (Shakura & Sunyaev 1973). Bondi accretion (Bondi 1952), meanwhile, implies very low angular momentum accretion – low enough that the gas encounters no angular momentum barrier on its way to the BH. Realistically, spherically symmetric Bondi accretion never occurs in nature, but three-dimensional (3D) time-dependent numerical simulations reveal that Bondi accretion can occur along the rotation axis of hot mode accretion flows (Janiuk, Proga & Kurosawa 2008).

Accretion discs themselves are thought to form from ‘Bondi-like’ accretion, in which a flow with low angular momentum accretes approximately spherically within the Bondi radius until it circularizes upon encountering a centrifugal barrier. Proga & Begelman (2003a) first explored this numerically by solving a specific instance of the ‘inhomogeneous Bondi problem’, adopting simple prescriptions for the azimuthal velocity, v_ϕ(r = R_o, θ), i.e. for the θ-component of angular momentum (l_θ). Their simulations confirmed the intuitive analytic results of Abramowicz & Zurek (1981): accretion at the Bondi rate can only occur if the gas does not encounter a centrifugal barrier, while hot mode accretion with a greatly reduced accretion rate results otherwise (unless the effective viscosity is quite large; see Narayan & Fabian 2011). Subsequent works have considered alternative prescriptions for angular momentum (e.g. Krumholz, McKee & Klein 2005), magnetohydrodynamic (MHD) evolution (Igumenshchev, Narayan & Abramowicz 2003; Proga & Begelman 2003b), and general relativistic effects for spinning BHs (Palit, Janiuk & Sukova 2019). In the last decade there have been numerous studies of low angular momentum accretion that include radiative feedback processes (e.g. Ciotti et al. 2009, 2010,2017; Kurosawa & Proga 2009; Park & Ricotti 2011; Li, Ostriker & Sunyaev 2013; Yuan et al. 2018; Bu & Yang 2019; Gan et al. 2019; Yoon et al. 2019), and this Bondi-like modelling framework has shown promise for reproducing the observed X-ray luminosity across many orders of magnitude in BH mass in early-type galaxies (e.g. Li et al. 2018; Pellegrini et al. 2018).

This Letter points to a need to reconsider the notion that accretion discs themselves can form via Bondi-like accretion. In general, the inhomogeneous Bondi problem has boundary conditions (BCs) defined on a spherical surface with a non-zero surface gradient (or surface curl) of density or pressure in addition to velocity. In this Letter, we investigate the other ‘simple’ axisymmetric problem: we consider a density BC with |$\partial$|ρ/|$\partial$|θ ≠ 0 rather than breaking spherical symmetry by considering |$\partial$|v_ϕ/|$\partial$|θ ≠ 0 at R_o. A non-uniform density distribution at R_o gives rise to pressure gradients in the flow that result in a non-zero v_θ. The meridional flows that develop are not at all analogous to the azimuthal flows in simulations with small non-zero l_θ. The streamlines have the opposite curvature, meaning an inflowing streamline turns into an outflowing one rather than circling the BH. Inflow–outflow solutions mediated by discs have long been speculated to occur due to the fact that the Bernoulli function is positive (Blandford & Begelman 1999, 2004; Begelman 2012), but here we are simply dealing with internal angular momentum transport without a disc.

This Letter is organized as follows. In Section 2, we show explicitly that the inhomogeneous Bondi problem can be considered an angular momentum transport problem, and we derive the critical density gradient that will result in inflow–outflow solutions. In Section 3, we present our simulation results, and in Section 4, we discuss the overall flow properties qualitatively.

2 THE INHOMOGENEOUS BONDI PROBLEM

The role of angular momentum in this problem is perhaps best appreciated through a quote from Kulsrud’s textbook (Kulsrud 2005). Referring to inviscid flows, he states that ‘steady state accretion cannot occur unless there is a magnetic torque that removes angular momentum’. First defining the total angular momentum in a general volume V by |$\boldsymbol {H} = \int _V \rho \, \boldsymbol {\ell }\, \mathrm{ d}V$|⁠, where |$\boldsymbol {\ell }= \boldsymbol {r}\times \boldsymbol {v}$| is the specific angular momentum of a fluid element located a distance |$\boldsymbol {r}$| away from some origin O, Kulsrud derives an equation for |$\mathrm{ d}\boldsymbol {H}/\mathrm{ d}t$| starting from the momentum equation in an ideal magnetohydrodynamic (MHD):

$$\begin{eqnarray*} \frac{\partial \left(\rho \boldsymbol {v}\right)}{\partial t} + {\nabla } \cdot \left(\rho \boldsymbol {v} \boldsymbol {v} + p\, \mathbb {I} + \boldsymbol {T}_\mathrm{ B}\right) = \rho \boldsymbol {g}, \end{eqnarray*}$$

(1)

where p is the gas pressure, |$\boldsymbol {T}_\mathrm{ B} =B^2/8\pi \mathbb {I} - \boldsymbol {B}\boldsymbol {B}/4\pi$| is the magnetic stress tensor (with |$\mathbb {I}$| the unit tensor and |$\boldsymbol{B}$| the magnetic field), and |$\boldsymbol {g}$| is the gravitational force. Namely, substituting for |$\partial \rho \boldsymbol {v}/\partial t$| in the expression for |$\mathrm{ d}\boldsymbol {H}/\mathrm{ d}t$| using equation (1), Kulsrud arrives at

$$\begin{eqnarray*} \frac{\mathrm{ d}\boldsymbol {H}}{\mathrm{ d}t} = -\int _S \rho \, \boldsymbol {\ell }\, \boldsymbol {v}\cdot \mathrm{ d}\boldsymbol {S} + \int _S \boldsymbol {r}\times \boldsymbol {B} \frac{\boldsymbol {B}}{4\pi } \cdot \mathrm{ d}\boldsymbol {S}. \end{eqnarray*}$$

(2)

Here, S is a spherical surface enclosing V. The only assumption used to derive equation (2) is that gravity is a central force, |$\boldsymbol {g} = -\nabla \Phi$|⁠. Hence, we see that in the absence of any viscosity, steady state accretion (⁠|$\mathrm{ d}\boldsymbol {H}/\mathrm{ d}t = 0$|⁠) can only be achieved if the above integrals are equal. This is the logic underlying the above quote. However, spherical, unmagnetized accretion (i.e. the Bondi solution) is the very special case |$\boldsymbol {B} = 0$| and |$\boldsymbol {\ell } = 0$|⁠, meaning the equality of equation (2) is trivially satisfied. The only exception to Kulsrud’s statement involves unmagnetized cases in which the integral |$\int _S \rho \, \boldsymbol {\ell }\, \boldsymbol {v}\cdot \mathrm{ d}\boldsymbol {S}$| vanishes due to internal dynamics in flows with non-zero specific angular momentum. Such is the nature of steady state solutions to the inhomogeneous Bondi problem: every spherical surface in V must satisfy |$\int _S \rho \, \boldsymbol {\ell }\, \boldsymbol {v}\cdot \mathrm{ d}\boldsymbol {S} =0$|⁠.

2.1 The critical density gradient

A BC with |$\partial$|ρ/|$\partial$|θ ≠ 0 at R_o leads to l_ϕ ≠ 0, and we therefore seek steady state solutions with v_ϕ = 0 but v_θ ≠ 0 at R_o. To identify a consistent set of BCs, in Appendix A we derive the basic equations governing steady state solutions. In particular, equation (A4) governs arbitrary axisymmetric BCs, and for a homentropic BC (⁠|$\partial$|s/|$\partial$|θ = 0) with v_ϕ = 0 it reads

$$\begin{eqnarray*} l_\phi = -R_\mathrm{ o}\frac{\partial h/\partial \theta + R_\mathrm{ o}\, v_\mathrm{ r} \, \partial v_\theta /\partial r }{v_\mathrm{ r} + \partial v_\theta /\partial \theta } . \end{eqnarray*}$$

(3)

For solutions obeying a standard outflow BC |$\partial$|v_θ/|$\partial$|r = 0, we can place a bound on the maximum density gradient that will result in pure inflow.¹ To accrete from large distances directly on to the BH, gas must have a specific angular momentum |$|\boldsymbol {\ell }| \lesssim 2 R_\mathrm{ s}\, c$| (Abramowicz & Zurek 1981), where R_s = 2GM_bh/c² is the Schwarzschild radius for a BH with mass M_bh. Enforcing this bound using equation (3), we find

$$\begin{eqnarray*} \left| \frac{\partial \ln \rho _{\rm {bdy}}}{\partial \theta } \right| \lesssim \frac{4R_\mathrm{ B}}{R_\mathrm{ o}} \frac{T_0}{T(R_\mathrm{ o},\theta)} \left| \frac{v_\mathrm{ r} + \partial v_\theta /\partial \theta }{c} \right|_{r = R_o} , \end{eqnarray*}$$

(4)

where |$R_\mathrm{ B} = GM_{\rm {bh}}/c_{\mathrm{ s},0}^2$| is the Bondi radius (with c_{s, 0} = (γp₀/ρ₀)^1/2 the adiabatic sound speed for a reference temperature |$T_0 = k_\mathrm{ B}^{-1} \bar{m} p_0/\rho _0$| at R₀). We have arrived at our main result: the quantity on the right-hand side is in general a tiny number. In the limiting case R_o → ∞, v_r(R_o) → 0, and |$\partial$|v_θ/|$\partial$|θ = 0 considered by Bondi (1952), ||$\partial$| ln ρ_bdy/|$\partial$|θ| can be infinitesimally small and yet still qualitatively change the nature of the Bondi solution by turning part of the inflow into an outflow, as we now demonstrate.

3 SIMULATIONS

Using athena++(version 19.0; Stone et al. in preparation), we solve the equations of adiabatic hydrodynamics, i.e. the force equation given in equation (1) with |$\boldsymbol {B} = 0$|⁠, the continuity equation |$\partial \rho /\partial t + \nabla \cdot \rho \boldsymbol {v} = 0$|⁠, and the entropy equation |$\partial s/\partial t + \boldsymbol {v}\cdot \nabla s = 0$|⁠, where s = c_v  ln p/ρ^γ is the specific entropy (with c_v the specific heat at constant volume and γ the adiabatic index). We examine solutions with γ = 5/3, implementing the Paczyńsky & Wiita (1980) pseudo-Newtonian potential, ϕ = −GM_bh/(r − R_s), to place the sonic point outside R_s. We perform 2D axisymmetric calculations in spherical polar coordinates (r, θ, ϕ) on a uniformly spaced polar grid with N_θ = 270 zones from 0 to |$\pi$|/2 and a logarithmic radial grid with N_r = 400 zones and dr_{i + 1}/dr_i = 1.019 from r_in = 10R_s to r_out = R_o. We apply reflecting BCs at θ = 0, |$\pi$|/2 and ‘constant gradient BCs’ (⁠|$\partial$|q/|$\partial$|r = constant) in all primitive variables at both the inner and outer radius unless otherwise noted. Finally, we use the third-order time integration scheme and the Harten–Lax–van Leer–Contact (HLLC) Riemann solver.

Solutions to the classic Bondi problem depend only on γ and R_B, whereas solutions with a boundary at finite R_o also depend on the value of R_o/R_B. This dependence has been recently explored by Samadi, Zanganeh & Abbassi (2019) and the solution procedure has been detailed by Waters & Proga (2012) in the context of Parker winds. Our initial conditions (ICs) are v_θ = v_ϕ = 0, while v_r, ρ, and p are the semi-analytic solutions to this generalized 1D Bondi problem under a BC ρ_bdy(R₀, θ) given below; the code for these ICs is publicly available at github.com/trwaters/bondiparker. We set R_o = 8R_B, noting that the outflow rates we find will have a dependence on this parameter since gas becomes less bound at larger radii (see Section 4).

We allow for an arbitrary density distribution to be specified along R_o using the function

$$\begin{eqnarray*} \rho _{\rm {bdy}}(\theta) = \rho _0\, \mathcal {N} f(\theta) , \end{eqnarray*}$$

(5)

where f(θ) is a distribution function and |$\mathcal {N}^{-1} = (1/2) \int _0^{\pi} f(\theta) \sin \theta \, \mathrm{ d}\theta$| is a normalization integral. The function f is normalized to recover the classic Bondi solution when f(θ) = 1, i.e. |$\int _0^{\pi} f(\theta)\, \mathrm{ d}\theta = \pi$|⁠. The normalization integral is chosen to satisfy |$M_\mathrm{ R} = \int \rho _{\rm {bdy}}(\theta)\, \mathrm{ d}V$| (⁠|${\approx} 2\pi \rho _0 \mathcal {N} \Delta r R_\mathrm{ o}^2 \int _0^{\pi} f(\theta) \sin \theta\, \mathrm{ d}\theta$| for Δr ≪ R_o), where |$M_\mathrm{ R} = 4\pi R_\mathrm{ o}^2 \Delta r \, \rho _0$| is the mass of a gas reservoir (a shell of matter extending from R_o to R_o + Δr that is constantly replenished) that is implicitly assumed to exist when applying this BC. By requiring that the mass of the reservoir is the same across all density distributions, we can properly compare the resulting accretion and outflow rates with the Bondi rate.

We consider the bell-shaped function f(θ) = (1 − δ cos ²θ)/(1 − δ/2), where δ is the additional free parameter governing our 2D solutions that controls the magnitude of the density inhomogeneity. For this choice, we can evaluate the left-hand side of equation (4) after averaging this bound over [0, |$\pi$|/2] to find |$(2/\pi) \int _0^{\pi /2}|\partial \ln f(\theta)/\partial \theta |\, \mathrm{ d}\theta = (2/\pi)|\ln (1 - \delta)| \approx 2\delta /\pi$| (provided δ ≪ 1); equation (4) becomes

$$\begin{eqnarray*} \delta \lesssim \frac{2\pi }{R_\mathrm{ o}/R_\mathrm{ B}} \left\langle \left| \frac{v_\mathrm{ r}/c + \partial (v_\theta /c)/\partial \theta }{T/T_0} \right| \right\rangle _{r = R_\mathrm{ o}}, \end{eqnarray*}$$

(6)

where 〈 · 〉 denotes an average from [0, |$\pi$|/2].

We present results for five runs: three with a constant entropy BC (A2, A3, and A4) and two with a constant pressure BC (B2 and B3; the number denotes the density contrast, −log δ). These BCs are implemented by assigning different pressures as a function of θ at R_o; we use p(R_o, θ) = p₀[ρ_bdy(θ)/ρ₀]^γ for runs A2, A3, and A4 and p(R_o, θ) = p₀ for runs B2 and B3, where p₀ = GM_bhρ₀/γR_B.

3.1 Solutions with constant entropy BCs

Table 1 summarizes the properties of our five runs. For solutions with homentropic BCs (runs A2, A3, and A4), the bound given by equation (6) is first formally satisfied for δ slightly below 10⁻⁴, meaning δ/δ_c falls below unity, δ_c being the value of the right-hand side of equation (6). Note that δ_c must be evaluated a posteriori since v_r(R_o, θ), |$\partial$|v_θ(R_o, θ)/|$\partial$|θ, and T(R_o, θ) are determined as part of a self-consistent solution. Table 1 shows that run A4 has a very weak outflow, consistent with δ/δ_c being somewhat larger than 1.

Table 1.

Properties of steady state solutions based on output at time 250R_B/c_{s, 0}. Runs A2, A3, and A4 (B2 and B3) are for constant specific entropy (pressure) along θ at R_o. The value of δ_c is found by evaluating the right-hand side of equation (6). The rates |$\dot{M}_{\rm {in}}$| and |$\dot{M}_{\rm {out}}$| are measured at R_o in units of |$\dot{M}_\mathrm{ B}$| and the accretion rate is |$\dot{M}_{\rm {A}} = \dot{M}_{\rm {in}} - \dot{M}_{\rm {out}}$|⁠. Multiplying |$\overline{M}_\mathrm{ o}$|⁠, the mean outflow Mach number at R_o (where the overbar denotes averaging only over outflowing zones), by (R_o/2R_B)^1/2 = 2 converts to velocity in units of the escape speed at R_o. The gas is unbound since the Bernoulli parameter at R_o is |$\mathrm{ Be} = 1.37c_{\mathrm{ s},0}^2 > 0$| for all of these runs; Be is only a weak function of radius, obeying |$\partial$|Be/|$\partial$|r = ℓ_ϕω_ϕ/r (see equation A2), where ω_ϕ is the ϕ-component of vorticity.

γ = 5/3, RB = 103Rs, Ro = 8RB
Run	δ	δ/δc	|$\dot{M}_{\rm {in}}$|	|$\dot{M}_{\rm {out}}$|	|$\overline{M}_\mathrm{ o}$|
A2	10−2	7.8	7.15	6.22	0.14
A3	10−3	3.6	2.22	1.26	0.032
A4	10−4	1.7	0.98	8.3 × 10−4	7.5 × 10−4
B2	10−2	–	1.34	0.36	0.013
B3	10−3	–	0.98	0.0	–

Steady state solutions for runs A2 and A3 are plotted in Fig. 1. The result of violating the angular momentum bound as measured by δ/δ_c is an increasingly strong outflow. The top panels compare the inflow, outflow, and net accretion rates |$(\dot{M}_{\rm {in}},\dot{M}_{\rm {out}},\dot{M}_{\rm {A}})$| as a function of radius. At R_o, run A3 has |$\dot{M}_{\rm {out}} \approx 1.3\dot{M}_{\rm {B}}$|⁠, while A2 has |$\dot{M}_{\rm {out}} \approx 6\dot{M}_{\rm {B}}$| (see the red lines). The net accretion rate (dashed black line) remains constant at nearly |$\dot{M}_{\rm {B}}$|⁠, and this is possible in a steady state only if |$\dot{M}_{\rm {in}} - \dot{M}_{\rm {out}} = \dot{M}_{\rm {A}} \approx \dot{M}_{\rm {B}}$| at every radius.

Figure 1.

Steady state dynamics of our 2D solutions. Top panels: inflow, outflow, and net accretion rates as a function of radius. Middle panels: log-scale colour maps of the ratio of the lateral pressure force and gravity, |r⁻¹|$\partial$|p/|$\partial$|θ/(GMρ/r²)|, with the velocity field overplotted. Red curves are contours of v_r = 0 and black curves are contours of ℓ_ϕ = 2R_sc. The inner border along these curves approximates the true ‘Bondi surface’ since regions interior to both curves have purely inflowing streamlines. This is made clear by plotting closely spaced streamlines around the inflow/outflow turning point (white curves in runs A2, A3, and B2). Bottom panels: log-scale colour maps of angular momentum (ℓ_ϕ) in units of R_sc zoomed-in on the region within the Bondi radius (⁠|$R_\mathrm{ B} = 10^3\, R_\mathrm{ s}$|⁠) as marked by the rectangles above. The velocity field is again overplotted. Note that ℓ_ϕ is less than R_sc and yet is a strongly varying function of (r, θ) despite other flow fields being closely spherically symmetric. Animations of these runs are viewable at https://trwaters.github.io/OutflowsFromInflows/.

The only force acting in the θ-direction that can drive outflow by turning streamlines around is the gas pressure force −r⁻¹|$\partial$|p/|$\partial$|θ. The middle panels in Fig. 1 show maps of the ratio of this force to gravity. Notice it is significant mainly in the outflowing regions at large radii. It changes sign between 4R_B and 6R_B, as it must transition from turning the flow outward to inward. The red and black contours denote where v_r = 0 and ℓ_ϕ = 2R_sc, respectively. Together these contours define an effective ‘Bondi surface’: the region internal to both contours has purely inflowing streamlines, while streamlines are outflowing in the external regions. This is seen clearly by identifying streamlines where the flow ‘turns over’ (white curves); this turnover point closely follows the red and black contours as δ is changed.

The bottom panels are maps of ℓ_ϕ in units of R_sc, zoomed-in on the region within 1R_B where the accretion is nearly spherical. These plots reveal that the flow is stratified in angular momentum, with ℓ_ϕ increasing outward. This dynamics is somewhat surprising, as gas at r ≲ R_B/2 has ℓ_ϕ ≪ R_sc, so intuitively angular momentum need not be transported outward for it to accrete. All that is required of these solutions from Section 2 is that they satisfy |$\int _S \rho \, \boldsymbol {\ell }\, \boldsymbol {v}\cdot \mathrm{ d}\boldsymbol {S} =0$|⁠, and since ρ and v_r are nearly spherically symmetric at r < R_B, this integral is approximately |$2\pi \rho v_\mathrm{ r} r^2 \int _0^{\pi} \ell _\phi \sin \theta \, \mathrm{ d}\theta$|⁠. The trivial way to obtain |$\int _0^{\pi} \ell _\phi \sin \theta \, \mathrm{ d}\theta = 0$| for ℓ_ϕ ≠ 0 is for ℓ_ϕ to be an odd function in the range [0, |$\pi$|], and this is automatic with an axis of symmetry at θ = |$\pi$|/2.

3.2 Solutions with constant pressure BCs

Constant pressure BCs are arguably more physical than constant entropy BCs, but in this case we can no longer derive a bound on |$\partial$| ln ρ_bdy/|$\partial$|θ. This is because equation (A4) now gives

$$\begin{eqnarray*} l_\phi = -R_\mathrm{ o}\frac{\partial h/\partial \theta - T\partial s/\partial \theta + R_\mathrm{ o}\, v_\mathrm{ r} \, \partial v_\theta /\partial r }{v_\mathrm{ r} + \partial v_\theta /\partial \theta } , \end{eqnarray*}$$

(7)

and the first two terms in the numerator are simply ρ⁻¹|$\partial$|p/|$\partial$|θ by the basic thermodynamic relation |$\rho ^{-1}\mathrm{ d}p = \mathrm{ d}h - T\, \mathrm{ d}s$|⁠. Thus, if we choose to assign a constant pressure at R_o, we are left with only an implicit relationship between ℓ_ϕ and |$\partial$|ρ/|$\partial$|θ through |$\partial$|v_θ/|$\partial$|θ via the continuity equation. Moreover, to allow for the possibility that ℓ_ϕ ≠ 0 at R_o when |$\partial$|p/|$\partial$|θ = 0, we apply a constant gradient BC to v_θ instead of the zero gradient BC (⁠|$\partial$|v_θ/|$\partial$|r = 0) used in Section 3.1.

Because a lateral gas pressure force (⁠|$=-R_\mathrm{ o}^{-1}\partial p/\partial \theta$|⁠) accompanies a constant entropy BC but not a constant pressure BC at R_o, we expect the outflow to be weaker for the same density BC. This is indeed the case, as Table 1 and the right two sets of panels in Fig. 1 show. Notice there is no black contour for runs B2 and B3 because all of the angular momentum in the domain is less than 2R_sc. Nevertheless, as a consequence of the density gradient |$\partial$|ρ/|$\partial$|θ at R_o, a lateral pressure force still builds up enough to cause a weak outflow in run B2 (qualitatively similar to that of run A3). We verified that an outflow as strong as in run A2 under constant pressure BCs can be obtained by increasing δ.

4 DISCUSSION

The adiabatic solutions we find overall consist of an inner purely inflowing region with ℓ_ϕ ≲ R_sc that is engulfed in a much larger region containing all of the (still low angular momentum) flow with ℓ_ϕ > R_sc. Rather than circularize, as would happen in solutions where the angular momentum is due to v_ϕ, the flow with ℓ_ϕ > R_sc forms an inflow–outflow region outside an effective ‘Bondi surface’, this being how angular momentum is transported outward. In recent papers following up an investigation by Hernández et al. (2014) using imposed density gradients as a mechanism for launching hydrodynamic jets, Aguayo-Ortiz, Tejeda & Hernandez (2019) and Tejeda, Aguayo-Ortiz & Hernandez (2019) found an incompressible analytic solution that captures the inflow–outflow topology of adiabatic solutions. However, these authors do not discuss the role of angular momentum, which is essential to understanding our results.

Our simulations show that the accretion rate on to the BH remains constant at nearly |$\dot{M}_\mathrm{ B}$|⁠, in agreement with the findings of Tejeda et al. (2019), although we do not expect this result to hold for simulations with a non-zero v_ϕ at R_o. We therefore caution against concluding that the Bondi formula provides an adequate estimate for the accretion rate on to ‘BH particles’ in cosmological simulations. More relevant for the construction of improved subgrid models of BH growth is our result that the outflow rate will always exceed |$\dot{M}_\mathrm{ B}$| for realistic density contrasts; this finding is in agreement with the 3D simulations of stellar wind fed accretion on to Sgr A* by Ressler et al. (2018). We note that spectral energy distribution (SED) fitting results for Sgr A* favour models in which the outflow transitions from weak to strong as radius increases (Ma et al. 2019), consistent with the dynamics we find here.

For a mean outflow velocity given by |$\overline{M}_\mathrm{ o} c_{\mathrm{ s},0}$| (with |$\overline{M}_\mathrm{ o}$| the mean Mach number of the outflow at R_o), the outflow rate should scale as |$\dot{M}_{\rm {out}} = \dot{M}_\mathrm{ B} \overline{M}_\mathrm{ o} (R_\mathrm{ o}/R_\mathrm{ B})^2$|⁠. In addition to the runs presented, which have R_o/R_B = 8, we examined variants of run A2 with R_o/R_B halved and doubled to verify this scaling. It is important to confirm that this scaling holds in Bondi–Hoyle-type simulations with an imposed wind (the properties of which determine the conditions along R_o), as well as to determine if the associated feedback can in turn limit the accretion rate by affecting the mass reservoir. We note that for supersonic motion through an inhomogeneous cloud, as considered in a few studies of 3D Bondi–Hoyle accretion (e.g. Ruffert 1999; Xu & Stone 2019), the imposed ram pressure may be too high and/or the density gradients too low to allow inflow–outflow regions to develop.² For marginally supersonic motion, however, this outflow dynamics should occur, and it is unclear what effect it will have on the instabilities accompanying Bondi–Hoyle accretion.

ACKNOWLEDGEMENTS

TW thanks Ya-Ping Li and Josh Dolence for helpful discussions and Sean Ressler for private correspondence. TW and AA are partially supported by the LANL LDRD Exploratory Research Grant 20170317ER. DP acknowledges support for program number HST-AR-14579.001-A provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. HL is supported by NASA/ATP and a LANL LDRD.

Footnotes

1

The enthalpy of an ideal gas is |$h = c_\mathrm{ s}^2/(\gamma - 1) = \gamma p \rho ^{-\gamma } \rho ^{\gamma -1}/(\gamma - 1)$|⁠, with |$c_\mathrm{ s}^2 = \gamma p/\rho$|⁠. Thus, taking the partial derivative |$\partial$|h/|$\partial$|θ at constant entropy gives |$\partial h/\partial \theta = c_\mathrm{ s}^2 \partial \ln \rho /\partial \theta$|⁠.

2

Specifically, taking |$R_\mathrm{ o} = GM_{\rm {bh}}/v_{\rm {\infty }}^2$| (the Hoyle–Lyttleton radius for relative motion through a cloud at speed v_∞), the quantity |$(R_\mathrm{ o}/R_\mathrm{ B})^2 = 1/M_{\rm {in}}^4$|⁠, where M_in = v_∞/c_{s, 0} is the characteristic Mach number of the mass reservoir, suggesting that inflow–outflow regions with |$\dot{M}_{\rm {out}} > \dot{M}_\mathrm{ B}$| will develop only for M_in ≲ 1.

APPENDIX A: RADIAL AND TRANSVERSE FLOW EQUATIONS

For simulations performed in spherical coordinates, it is useful to derive equations governing the flow on and normal to any spherical surface in a steady state; these follow from Crocco’s theorem (see e.g. Rezzolla & Zanotti 2013):

$$\begin{eqnarray*} \nabla \mathrm{ Be} = \boldsymbol {v} \times \boldsymbol \omega + T\nabla s + \boldsymbol {f}_\mathrm{ b}, \end{eqnarray*}$$

(A1)

where Be ≡ h + v²/2 + Φ is the Bernoulli function, |$\boldsymbol \omega \equiv \nabla \times \boldsymbol {v}$| the vorticity, and |$\boldsymbol {f}_\mathrm{ b}$| the sum of both magnetic and radiation forces. Note that Crocco’s theorem governs both adiabatic and non-adiabatic solutions and that the generalized Bernoulli’s theorem follows by dotting equation (A1) with |$\boldsymbol {v}$| to give |$\nabla \mathrm{ Be} =T\nabla s + \boldsymbol {f}_\mathrm{ b}$| along streamlines. By instead dotting equation (A1) with |$\boldsymbol {r} = r\hat{r}$|⁠, we obtain the radial flow equation

$$\begin{eqnarray*} \boldsymbol {r} \cdot \nabla \mathrm{ Be} = \boldsymbol {\ell } \cdot \boldsymbol \omega + T\boldsymbol {r}\cdot \nabla s + \boldsymbol {r}\cdot \boldsymbol {f}_\mathrm{ b}. \end{eqnarray*}$$

(A2)

In the Bondi solution, all terms on the right-hand side are zero and |$\boldsymbol {v}\parallel \boldsymbol {r}$|⁠, so this reduces to the classical Bernoulli theorem, namely Be = constant along streamlines.

Taking the cross product of |$\boldsymbol {r}$| and equation (A1) gives

$$\begin{eqnarray*} \boldsymbol {r}\times \nabla \left(h + \frac{v^2}{2} \right) &=& (\boldsymbol {v}\cdot \boldsymbol \omega)\boldsymbol {r} - r\, v_\mathrm{ r}\boldsymbol \omega - \boldsymbol \omega \times \boldsymbol {\ell } + T\boldsymbol {r}\nonumber\\ &&\times \,\nabla s + \boldsymbol {\tau }, \end{eqnarray*}$$

(A3)

where the identities |$\boldsymbol {r} \times (\boldsymbol {v} \times \boldsymbol \omega) = \boldsymbol {v}\times (\boldsymbol {r}\times \boldsymbol \omega) - \boldsymbol \omega \times \boldsymbol {\ell }$| and |$\boldsymbol {v}\times (\boldsymbol {r}\times \boldsymbol \omega) = (\boldsymbol {v}\cdot \boldsymbol \omega)\boldsymbol {r} - (\boldsymbol {v}\cdot \boldsymbol {r})\boldsymbol \omega$| were used. Here, |$\boldsymbol {\tau } = \boldsymbol {r}\times \boldsymbol {f}_\mathrm{ b}$| is the torque from any magnetic or radiation forces. Note that the radial component of this equation is an identity. In axisymmetry, the θ-component is τ_θ/r = v_rω_θ − v_θω_r, and both ω_r and ω_θ are zero when v_ϕ = 0, leading to an inconsistency if τ_θ ≠ 0. This is the intuitive statement that steady state solutions with v_ϕ = 0 are impossible if there is any net applied torque in the poloidal plane, for meridional circulation will just continually build up. From the ϕ-component, we can arrive at an expression for ℓ_ϕ = rv_θ for axisymmetric flows that holds on every spherical surface:

$$\begin{eqnarray*} l_\phi = r\frac{v_\phi ^2\cos \theta - r v_\mathrm{ r} \sin \theta \, \partial v_\theta /\partial r + S(\boldsymbol {r})\sin \theta }{(v_\mathrm{ r} + \partial v_\theta /\partial \theta)\sin \theta }, \end{eqnarray*}$$

(A4)

where |$S(\boldsymbol {r}) = T\partial s/\partial \theta - \partial h/\partial \theta + \tau _\phi$|⁠. Notice from equation (A4) that stagnation points (where |$\boldsymbol {v} = 0$|⁠) can only occur at locations |$\boldsymbol {r}_\mathrm{ s}$| satisfying |$S(\boldsymbol {r}_\mathrm{ s}) = 0$|⁠, meaning there is a lateral force balance, |$\partial$|p/|$\partial$|θ = ρτ_ϕ, while there is simultaneously a radial force balance by equation (A2).