Simulation of a compact object with outflows moving through a gaseous background

ABSTRACT

A compact object moving relative to surrounding gas accretes material and perturbs the density of gas in its vicinity. In the classical picture of Bondi–Hoyle–Lyttleton accretion, the perturbation takes the form of an overdense wake behind the object, which exerts a dynamical friction drag. We use hydrodynamic simulations to investigate how the accretion rate and strength of dynamical friction are modified by the presence of outflow from the compact object. We show that the destruction of the wake by an outflow reduces dynamical friction, and reverses its sign when the outflow is strong enough, in good quantitative agreement with analytic calculations. For a strong isotropic outflow, the outcome on scales that we have simulated is a negative dynamical friction, i.e. net acceleration. For jet-like outflows driven by reprocessed accretion, both the rate of accretion and the magnitude of dynamical friction drop for more powerful jets. The accretion rate is strongly intermittent when the jet points to the same direction as the motion of the compact object. The dynamical effects of outflows may be important for the evolution of compact objects during the common envelope phase of binary systems, and for accreting compact objects and massive stars encountering active galactic nucleus discs.

accretion, accretion discs, stars: neutron, stars: winds, outflows

1 INTRODUCTION

When a massive object moves through a cloud of gas, its gravity pulls ambient gas towards it leading to accretion. In the standard picture of Bondi–Hoyle–Lyttleton (BHL) accretion (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Edgar 2004), the flow is axisymmetric and gas whose unperturbed streamlines pass inside the Bondi radius ¹

$$\begin{eqnarray*} R_{ B} = \frac{2GM}{c_s^2+V_*^2} \end{eqnarray*}$$

(1)

is accreted on to the object. Here, M and V_* are the mass and velocity of the object and c_s is the sound speed of the ambient gas. The total accretion rate, valid for both subsonic and supersonic cases, is given by Bondi (1952) as

$$\begin{eqnarray*} \dot{M}_B=\alpha \pi R_{ B}^2\rho _0 V_*\approx \frac{4\pi \alpha (GM)^2\rho _0}{(V_*^2+c_s^2)^{3/2}} \end{eqnarray*}$$

(2)

where ρ₀ is the ambient gas density. α is a constant of order unity first introduced by Bondi & Hoyle (1944) to account for the non-uniformity of gas velocity inside the Bondi radius due to the gravitational circulation. Bondi (1952) uses α = 1/2, but numerical calculations by Shima et al. (1985) suggest α = 1. Due to the long-range nature of gravity, gas residing outside the Bondi radius also experiences gravitational forces. As a consequence, incoming gas with an impact parameter larger than the Bondi radius is not accreted on to the star, but accumulates behind the star, forming a overdense wake. Gravity from the overdense wake exerts a drag opposite to the object’s direction of motion known as dynamical friction (Chandrasekhar 1943; Ostriker 1999). For the supersonic case, the acceleration of dynamical friction is given by

$$\begin{eqnarray*} a_{DF}=\frac{4\pi G^2 M\rho _0}{V_*^2}\ln \left[ \Lambda \left(1-\frac{1}{\mathcal {M}^2}\right)^{0.5}\right] \end{eqnarray*}$$

(3)

where |$\mathcal {M}=V_*/c_s$| is the Mach number of the object, and ln Λ = ln (b_max/b_min) is the Coulomb logarithm with b_max and b_min being the maximum and minimum impact parameters, respectively. ln Λ quantifies the spatial range over which gravity is important.² Numerical simulations are in generally good agreement with this simple picture (Blondin & Raymer 2012).

The standard analytic model of BHL accretion does not include any outflow from the central object. Outflows, however, are expected to be present in many of the physical situations where BHL accretion occurs. Massive stars can drive winds and previous studies (Soker 2016; Murray, Goyal & Chang 2018; Shiber et al. 2019; Hillel, Schreier & Soker 2019) have found that outflow may play an important role in the common envelope evolution of binary stars. In particular, powerful outflows are observed to accompany super-Eddington accretion on to compact objects in sources such as SS433 (Fabrika 2004). The influence of such outflows on accretion and dynamical friction remains to be fully explored. Gruzinov, Levin & Matzner (2019) considered the case where a strong isotropic wind is driven from the compact object. Their analytical calculations showed that when the strong wind shuts off the accretion, a bow shock is launched around the object. Instead of an overdense wake, an underdense hole is formed behind the bow shock. The gravitational acceleration of the object from the ambient gas becomes ‘negative dynamical friction’³ if the wind velocity V_w and the object’s velocity V_* satisfy u ≡ V_*/V_w < 1.71. In addition to this purely mechanical effect, thermal and radiative feedback can also be important (Taam, Fu & Fryxell 1991; Blondin et al. 1990; Kley, Shankar & Burkert 1995; Park & Ricotti 2013; Souza Lima et al. 2017). Park & Bogdanović (2017) studied the radiative feedback of accreting massive black holes and found that the overdense wake can be destroyed, reducing the dynamical friction and in some cases reversing its direction. Masset & Velasco Romero (2017) showed that diffusion of heat through an opaque medium surrounding the accreting object similarly leads to a heating force of opposite sign (and in some cases larger magnitude) than the dynamical friction. Related physical effects may be important for the migration of planets (and generically for embedded gravitating objects) in Keplerian disc flows (Masset 2017).

In this paper, we carry out a numerical study of a compact object with an outflow moving through uniform ambient gas. We focus on the case where the compact object is moving supersonically through the gas and study the effects of outflow on the accretion rate and dynamical friction in the case of both a beamed jet as well as an isotropic wind. To make the problem numerically tractable, we parametrize the properties of the outflow and do not self-consistently model the launching mechanism. The plan of the paper is as follows. In Section 2, we outline our numerical methods and set up for the simulations. In Section 3, we present results for the case of an isotropic wind and compare it to the analytical calculations of Gruzinov et al. (2019). The case of a beamed jet is presented in Section 4. The final section is devoted to discussions.

2 NUMERICAL SCHEMES

We use the ramses code (Teyssier 2002) to simulate a compact object moving in a 3D Cartesian box. ramses is a mature adaptive mesh refinement (AMR) code that includes self-gravity, sink particles (Dubois et al. 2010; Bleuler & Teyssier 2014), and radiation (Rosdahl et al. 2013). Previously, one of the authors and his collaborator used ramses to study star formation in a turbulent gas both without feedback (Murray et al. 2017) and with jet feedback (Murray et al. 2018). Here, we use the modifications for sink particles and jet and wind feedback for ramses that is described in Murray et al. (2018).

The sink particle module that we use in ramses is the same as described in Murray et al. (2017), but with a few crucial changes. First, we create the sink particle at the beginning of the simulation, not based on a Jeans criterion as in Murray et al. (2017). Second, we accrete all the gas above a (low) threshold density within two cells of the sink particle. ⁴ Similarly, we use the same jet and wind feedback module as Murray et al. (2018). For completeness, we briefly describe the jet and wind feedback module here. We inject mass, momentum and energy in a region that is 4–8 cells from the sink particle at the finest level. For the isotropic wind, we use a flat distribution in angles. For the jet case, we assume a Gaussian jet. Here, the injection region consists of a bi-cone with a user-defined opening angle of 90° (45° half-opening angle) about a chosen jet axis. We note that while the initial angle can be large, recollimation produces a narrower jet on large scales.

The 3D simulation box is set up in the rest frame of the compact object, where the compact object with mass M is placed at the centre of the box, as in MacLeod et al. (2017). In our simulations, the ambient ideal gas initially has constant density ρ₀ flowing along the +x direction with constant velocity V_* (the compact object is moving along −x direction in the lab frame). The gas obeys an ideal equation of state, with adiabatic index γ = 5/3. The boundaries of the computational box are set to be periodic in the y- and z-directions and ingoing/outgoing for the left/right boundary in the x-direction. We do not account for deflection of streamlines outside the box, and the box size and resolution are chosen so that L_B ≫ R_B to reduce the effects of gravity on the fluid near the boundary. Our grid spacing Δx is chosen to be smaller than R_B to resolve the relevant scale of the accretion process as well as the structure formed by the outflow (e.g. bow shock). In our typical runs, we use 256 cells for each side with one level of refinement around the compact object, for an effective resolution of 512³. We have also performed convergence tests to make sure that our simulations produce consistent results independent of the box size and resolution.

The compact object is treated as a sink particle in the simulation which accretes ambient gas with its own gravity, which is smoothed over 4 cells on the finest level. The simulation includes both the gravity of the compact object and self-gravity of the gas. However, as we use periodic boundary conditions for our gravity solver, the gravity from the uniform gas background is neglected, i.e. the zeroth Fourier mode of the gas density distribution is ignored. This affects the self-consistency of our results on the largest scales, as we discuss in Section 4.3. In addition, the compact object is accelerated as it absorbs the momentum of the gas it accretes. (Any outflow emits no net momentum.) We update the object’s mass and momentum as it accretes, thereby tracking the components of its net acceleration; however, we do not update its position, which is fixed at the box centre.

The relevant parameters and observables in our simulations are the ambient gas density ρ₀ and sound speed c_s, the velocity V_*, mass M and accretion rate |$\dot{M}$| of the compact object, the velocity V_w, and mass-loss rate |$\dot{m}_w$| of the outflow. We will present our results in dimensionless parameters

$$\begin{eqnarray*} \mathcal {M}\equiv \frac{V_*}{c_s} \quad \textrm {and}\quad u\equiv \frac{V_*}{V_w}. \end{eqnarray*}$$

(4)

The acceleration will be normalized by a_DF evaluated with Λ = 64 being the ratio of the half-box size over the smoothing length of the gravity. All other quantities will be normalized by the corresponding values from the Bondi accretion (α = 1) and we define the Bondi time-scale as

$$\begin{eqnarray*} t_B\equiv \frac{2GM}{\left(c_s^2+V_*^2\right)^{3/2}}. \end{eqnarray*}$$

(5)

We run each simulation for a time much longer than t_B to allow the system to relax.

3 ISOTROPIC WIND

We first simulate the case where the wind blowing from the star is isotropic with constant mass-loss rate |$\dot{m}_w$| and velocity V_w. The star is taken to be a M = 1 M_⊙ compact object and the ambient gas initially has constant density ρ₀ = 10⁻⁶ g cm⁻³ with pressure p₀ corresponding to a sound speed of c_s0 = 38 km s⁻¹. The fiducial parameters used in the simulation are set as V_* = 600 km s⁻¹ and V_w = 3000 km s⁻¹, corresponding to Mach number of |$\mathcal {M}\approx 15.8$| and u = 0.2. The Bondi values and dynamical friction for this setup are

$$\begin{eqnarray*} R_B &=& 7.4\times 10^{10}\,\,\mathrm{cm},\nonumber \\ \dot{M}_B &=& 0.016 \,\, M_\odot \,\mathrm{yr}^{-1},\nonumber \\ t_B &=& 1230\,\,\mathrm{s},\nonumber \\ a_{DF} &=& 0.15\,\,\mathrm{cm\,s}^{-2}. \end{eqnarray*}$$

We set the box size to be L_B = 10¹³ cm, much larger than R_B, and adopt the mass-loss rate |$\dot{m}_w=1 M_\odot$| yr⁻¹. Here, |$\dot{m}_w$| is much larger than the Bondi accretion rate |$\dot{M}_B$| to test the analytical solution given by Gruzinov et al. (2019), which neglected the influence of the object’s gravity on the gas flow.

Fig. 1 shows the gas density slice in the central x − y plane at t = 630t_B . The colour shows the gas density and the blue streamline follows the gas velocity. The strong wind outflow blocks accretion and redirects the flow, and forms a high-density bow shock surrounding an underdense cavity around the compact object. At the surface of the bow shock, the gas density exhibits corrugated patterns indicating the presence of the Kelvin–Helmholtz instability.

Figure 1.

Slice of gas density in the central x − y plane with isotropic wind at t = 630t_B. The streamline shows the gas velocity and the white line is the analytical solution from Wilkin (1996).

The white line in Fig. 1 shows the the analytical solution given by Wilkin (1996), which neglects gravitational effects and models the wind-ambient collision as inelastic, resulting in formation of a thin shell. The white line marks the position of the shell in the frame of the compact object, given by

$$\begin{eqnarray*} R_s = R_0\frac{\sqrt{3(1-\theta \cot \theta)}}{\sin \theta }, \end{eqnarray*}$$

(6)

in spherical coordinates for which the central compact object is at the origin, and the axis θ = 0 is in the negative-x direction. Here, |$R_0\equiv \sqrt{\dot{m}_w V_w/4\pi \rho _0 V_*^2}$| is the characteristic standoff distance. The overdensity created by the bow shock follows the white line but lies further away from the compact object since the adiabatic ambient shock separates from the contact discontinuity.

Fig. 2 shows the density, velocity, and pressure profiles along the x-axis at the centre of the computational box. The dashed line indicates the shock position from the analytical solution by Wilkin (1996). The position of the shock marked by the velocity jump is close to the dashed line, and the density increases ahead of the shock as gas flow from upstream accumulates there. Moreover, the reverse shock is close to the central compact object, resulting in a large pressure which also pushes the overdensity outward. Behind the compact object, the wind blows away the gas causing the density to drop while the velocity remains constant.

Figure 2.

Gas profile along the x-axis at the centre. Upper panel: gas density, middle panel: velocity parallel to the x-axis, and bottom panel: pressure. The dashed line is the analytical shock position from Wilkin (1996).

Gruzinov et al. (2019) derive an analytical approximation for the gravitational acceleration on the central compact object, using the matter distribution predicted by Wilkin (1996). In Wilkin’s solution, the density has constant value ρ₀ outside the bow shock and drops as |$\rho _0 u^2R_0^2/R^2$| from the central compact object inside the bow shock. The surface density of the bow shock follows

$$\begin{eqnarray*} \sigma = \rho _0 R_0^2\frac{\left[2u(1-\cos \theta)+R_s^2\sin ^2\theta /R_0^2\right]^2}{2R_s\sin \theta \sqrt{(\theta -\sin \theta \cos \theta)^2+\left(R_s^2/R_0^2-1\right)\sin ^4\theta }}. \end{eqnarray*}$$

(7)

The total gravitational force is pointing along the x-direction due to symmetry and has contributions from both the gas and the bow shock

$$\begin{eqnarray*} F_{\rm gas} = - 2\pi GM\int _0^\infty r^2 \mathrm{d}r \int _0^\pi \sin \theta \mathrm{d}\theta \frac{\rho _0}{r^2}\cos \theta , \end{eqnarray*}$$

(8)

$$\begin{eqnarray*} F_{\rm shock} = -2\pi GM \int _0^\pi J_\theta \mathrm{d}\theta R_s\sin \theta \frac{\sigma }{r^2} \cos \theta , \end{eqnarray*}$$

(9)

where J_θ is the Jacobian

$$\begin{eqnarray*} J_\theta \equiv \sqrt{\left(\frac{\partial R_s\sin \theta }{\partial \theta }\right)^2 + \left(\frac{\partial R_s\cos \theta }{\partial \theta }\right)^2}. \end{eqnarray*}$$

(10)

The net force is

$$\begin{eqnarray*} F &=& -\pi GM\rho _0 \int _0^\pi \mathrm{d}\theta \cos \theta \sin \theta R_s \nonumber \\ &&\times \left[ \frac{3}{2}\left(1+\frac{2u(1-\cos \theta)}{R_s^2\sin ^2\theta /R_0^2} \right)^2 - 2\left(1+\frac{u^2}{R_s^2/R_0^2}\right) \right]. \end{eqnarray*}$$

(11)

For u ≪ 1, the net force asympototes to

$$\begin{eqnarray*} F = -8.18 GM\rho _0 R_0 \end{eqnarray*}$$

(12)

which points to the same direction as the object’s motion.

Equation (11) assumes the ambient gas follows the distribution predicted by Wilkin (1996) towards the spatial infinity. However, in the realistic environment and numerical simulations, the matter distribution deviates from Wilkin's (1996) prediction on sufficiently large scales. (We discuss this issue further in Section 4.3.) We compute the net force due to gas within a sphere of radius |$\tilde{R}$|⁠: this can be derived by replacing the upper limit in the radial integral with |$\tilde{R}$|⁠.

$$\begin{eqnarray*} \int _0^{\tilde{R}} r^2 \mathrm{d}r \frac{\rho _0}{r^2} &=& \int _0^{R_s} \mathrm{d}r \rho _0 \frac{ u^2 R_0^2}{r^2} + \int _{R_s}^{\tilde{R}} \mathrm{d}r \rho _0 \nonumber \\ &=& \rho _0\left(\tilde{R}-R_s-\frac{u^2 R_0^2}{R_s}\right), \end{eqnarray*}$$

(13)

for |$\theta \lt \tilde{\theta }$| where the bow shock is inside the sphere. The angle |$\tilde{\theta }$| is defined to be where the bow shock crosses the sphere |$R_s(\tilde{\theta })=\tilde{R}$|⁠. For the range of polar angle where the bow shock is outside the sphere

$$\begin{eqnarray*} \int _0^{\tilde{R}} r^2 \mathrm{d}r \frac{\rho _0}{r^2} = \int _0^{\tilde{R}} \mathrm{d}r \rho _0\frac{ u^2 R_0^2}{r^2} = -\rho _0\frac{u^2 R_0^2}{\tilde{R}}. \end{eqnarray*}$$

(14)

The upper limit in the force from the shock should be be replaced by |$\tilde{\theta }$| in the angular integral since the bow shock outside the sphere will not contribute to the gravity. Therefore, the expression for the net gravity of gas within a sphere of radius |$\tilde{R}$| is

$$\begin{eqnarray*} F(\tilde{R}) &=& -\pi GM\rho _0 \int _0^{\tilde{\theta }} \mathrm{d}\theta \cos \theta \sin \theta R_s \nonumber \\ &&\times \left[ \frac{3}{2}\left(1+\frac{2u(1-\cos \theta)}{R_s^2\sin ^2\theta /R_0^2} \right)^2 - 2\left(1+\frac{u^2}{R_s^2/R_0^2}\right) \right] \nonumber \\ &&-\, \pi GM\rho _0\left(\tilde{R}+\frac{u^2 R_0^2}{\tilde{R}^2}\right) \sin ^2\tilde{\theta } . \end{eqnarray*}$$

(15)

Equation (15) asymptotes to equation (11) for large |$\tilde{R}$| at the order of |$\mathcal {O}(\tilde{R}^{-1})$|⁠. Fig. 3 shows the relative error of equation (15) as a function of |$\tilde{R}$| for various u, and shows this R⁻¹ scaling. For u < 0.5 and |$\tilde{R}\gtrsim 5 R_0$|⁠, equation (15) agrees with equation (11) for better than |$10{{\ \rm per\ cent}}$|⁠. On this basis, we would expect the net gravitational acceleration obtained from the simulation to be close to the value given be equation (11) as long as our computational box is larger than ∼5R₀ for u < 0.5. For the fiducial parameters we use, we have L_B/2R₀ = 7.7, and we should expect convergence to better than |$10{{\ \rm per\ cent}}$|⁠.⁵

Figure 3.

Relative error of equation (15) as a function |$\tilde{R}$| for various u ≡ V_*/V_w.

Fig. 4 shows the evolution of the gravitational acceleration for the compact object in the x-direction from the ambient gas. After an initial stage of increase where the wind flows out to form a steady bow shock, the acceleration of the compact object settles to oscillate around a steady value. The dotted line is the average acceleration for t > 100t_B yr and the dashed line is the theoretical calculation from Gruzinov et al. (2019). We can see that the two values are close to each other. The value from Gruzinov et al. (2019) is larger since their calculation assumes the overdensity sits at the shock, but our simulation suggests that the overdensity is farther away. We have also compared the acceleration obtained with high-resolution simulations and larger box size, with entirely consistent results.

Figure 4.

Evolution for the magnitude of the gravitational acceleration of the compact object in the x-direction from the ambient gas. The acceleration points in the −x direction and has negative values. The dotted line is the average acceleration for t > 100t_B and the dashed line is the theoretical calculation from Gruzinov et al. (2019).

Fig. 5 shows the gravitational acceleration for different values of u (varying V_w while holding V_* constant) with all other parameters kept fixed. The ratio |$L_B/2R_0=17.24\sqrt{u}\gt 5$| for all the cases considered in the plots. Black dots are simulation results which is the average acceleration after a steady state is reached. Red dots are analytical calculations (equation 11) and blue dots are the asymptotic expression for u ≪ 1 (equation 12). We can see that the simulation results are all within |$\sim 10{{\ \rm per\ cent}}$| of the analytical results, and that the asymptotic formula for u ≪ 1 tends to overestimate the acceleration for non-zero u.

Figure 5.

The dependence of acceleration on the wind speed V_w. Black dots are simulation results. Red dots are analytical results (equation 11) and blue dots are the asymptotic expression for u ≪ 1 (equation 12).

4 ACCRETION-POWERED JET

The numerical results presented in Section 3 is for the purpose of testing the analytical solution, where we use a fixed mass-loss rate |$\dot{M}$| larger than the Bondi accretion rate. In this section, we study a more physical scenario where a compact object moves on the background of ambient gas with accretion powered outflow. The compact object still has mass of 1 M_⊙, but launches beamed jet rather than isotropic wind. We set the star velocity to be V_* = 300 km s⁻¹ corresponding to Mach number |$\mathcal {M}\approx 7.9$|⁠. The Bondi values and the predicted dynamical friction force for this setup are

$$\begin{eqnarray*} R_B &=& 3 \times 10^{11}\,\,\mathrm{cm},\nonumber \\ \dot{M}_B &=& 0.13 \,\, M_\odot \,\mathrm{yr}^{-1},\nonumber \\ t_B &=& 9735\,\,\mathrm{s},\nonumber \\ a_{DF} &=& 0.6\,\,\mathrm{cm\,s}^{-2}.\nonumber \end{eqnarray*}$$

Again, we set the box size to be 10¹³ cm, much larger than R_B.

4.1 No outflow

As a benchmark for comparison, we first run simulations with no outflow. The supersonic motion of the object leaves an axis-symmetric conic wake behind. Fig. 6 shows the slice of gas density in the x − y plane. The velocity streamline of ambient gas is refracted towards the object by the gravitational drag until it reaches the wake. The overdense gas inside the wake acts gravitationally to decelerate the object which is the standard dynamical friction. The values of mass accretion rate and dynamical acceleration are consistent with conventional analytical values given by equations (2) and  (3) as shown in Figs 8 and 9 for the case of no outflow.

Figure 6.

Slice of gas density in the x−y plane for the simulation with no outflow at t = 80t_B.

4.2 Beamed Jet

For accretion powered outflow, the mass outflow rate is expected to be a fraction of the mass accretion rate. Therefore, we set the mass-loss rate |$\dot{m}_w$| to be a constant fraction

$$\begin{eqnarray*} f\equiv \frac{\dot{m}_w}{\dot{M}} \end{eqnarray*}$$

(16)

of the numerical mass accretion rate |$\dot{M}$|⁠. As discussed in Section 2, |$\dot{M}$| is calculated by accreting all the gas within two cells from the sink particle at the finest level. The shape of the outflow has a Gaussian profile spreading over a polar angle θ = 45° from the beam centre. We vary the parameter u by changing the outflow velocity V_w, while keeping V_* constant and make sure that the size of the consequent bow shock remains smaller than the box size of L_B = 10¹³ cm.

Fig. 7 shows the slices of the simulation where the outflow is beamed around the z-axis (at right angles to the incoming gas) with u = 0.03 (V_w = 10⁴ km s⁻¹) and f = 0.3. The slices are plotted in three planes passing through the box centre. In the y − z slice on the top, a bipolar outflow is observed. A shock front with gas overdensity is formed ∼4R_B in the z-direction from the central object. Close to the compact object, there is a underdense gas bubble along the z-axis, but overdense gas along the y-axis. In the x − z slice at the centre, a bow shock is clearly seen as the beamed outflow encounters the ambient gas. The jet sweeps up the gas, creating an underdense tunnel which extends along the z-axis close to the compact object but is bent in the x-direction by the ambient gas far away from the object. Compared to Fig. 1, the underdensity inside the bow shock less prominent, for two reasons: the outflow is now beamed, so it sweeps a smaller region; and the mass outflow rate is much smaller. Moreover, there is an overdense region between the underdense tunnels swept out by the jet. The x − y slice on the bottom is similar to Fig. 6 without outflow, indicating gas not swept by the jet still accumulates behind the star due to gravitational drag.

Figure 7.

Slices of gas density for three planes passing through the computation box centre at t = 80t_B after the simulation starts. The outflow is beamed around the z-axis perpendicular to the velocity of the compact object with u = 0.03 (V_w = 10⁴ km s⁻¹) and f = 0.3.

Fig. 8 shows the comparison of accretion rate evolution between simulation with and without the outflow. The mass of compact object increases as it accretes matter from the gas, and therefore the accretion rate should also change as M² in the simulation. However, since the accretion rate is low (0.1 M_⊙ yr⁻¹) and the simulation time is short, the change of mass is negligible (∼10⁻³). Therefore, the accretion rate appears to be stationary after an initial phase of relaxation. The measured accretion rate for the case without outflow is |$\dot{M}_0=0.7$| |$\dot{M}_B$|⁠, consistent with the BHL accretion theory. As the jet is turned on, the outflow disrupts inflowing gas and reduces the accretion rate on to the object. The mass accretion rate keeps dropping with the increasing outflow velocity of the gas.

Figure 8.

Comparison of mass accretion rate evolution for different simulations. The beamed jet is launched along z-axis with f = 0.3 and various with various u (various V_w, V_* is kept constant). The dashed lines show the average acceleration for t > 30t_B.

Fig. 9 shows the net acceleration of the star along its direction of motion. The acceleration is dominated by the gravitational force from the ambient gas, but also includes contributions (⁠|$\sim 10{{\ \rm per\ cent}}$|⁠) from the accreted momentum. Without the outflow, the gas accumulated behind the star acts as the source of the standard dynamical friction, and its magnitude is close to our estimation of a_DF. As the outflow velocity V_w increases, the dynamical friction is reduced by the underdense tunnels created by the jet. However, no stable negative dynamical friction is observed in our simulations. Only a transient stage of negative acceleration is observed during the early epoch for the case with strongest outflow u = 0.03 (V_w = 10⁴ km s⁻¹).

Figure 9.

Comparison of acceleration of the object in the x-direction for different simulations. The beamed jet is launched along z-axis with various u (various V_w, V_* is kept constant). The dashed lines show the average acceleration for t > 30t_B.

In addition to the scenario in which the jet launched perpendicular to the direction of the object’s motion, we have also simulated cases where the jet is parallel. This could be relevant if its orientation aligns with the spin of the compact object, which could point in any direction. Fig. 10 shows slices for the parallel case with u = 0.03 (V_w = 10⁴ km s⁻¹) and f = 0.3. There is a region of underdensity around the compact object seen in all three plots. Unlike the perpendicular case, there is no tunnel of underdensity observed.

Figure 10.

Slices of gas density for three planes passing through the computation box centre at t = 80t_B . The outflow is beamed around the x-axis perpendicular to the velocity of the compact object with u = 0.03 (V_w = 10⁴ km s⁻¹) and f = 0.3.

Fig. 11 shows the accretion rate evolution when the outflow is beamed in the x- direction (along the flow axis). Unlike a perpendicular jet, there is no steady state in this case: sporadic spikes in |$\dot{M}$| punctuate a low level of accretion. When the accretion rate is high, the mass-loss rate also increases which feeds back negatively on the accretion rate that in turn reduces the mass-loss rate. As the outflow suffers from the head-on collision with the ambient gas, the mixing of streams can excite gas instabilities and break the periodicity of the evolution leading to intermittent spikes. Signatures of the spikes in accretion rate are also present in the slice plots, where sporadic low-density patterns are ejected from the object which then propagate away in the overdense wake.

Figure 11.

Comparison of mass accretion rate evolution for different simulations. The beamed jet is launched along x-axis with f = 0.3 and various u (various V_w, V_* is kept constant). The dashed lines show the average acceleration for t > 30t_B.

The intermittent pattern is also present in the acceleration of the object as shown in Fig. 12. Even though the acceleration is more irregular, its values are still positive with no negative dynamical friction observed in the simulations.

Figure 12.

Comparison of acceleration of the object in the x-direction for different simulations. The beamed jet is launched along x-axis with various u (various V_w, V_* is kept constant). The dashed lines show the average acceleration for t > 30t_B.

Figs 13 and 14 show how the accretion rate and acceleration can change with the outflow efficiency f. With increasing f, both the mass accretion rate and the acceleration decrease whether the jet is along x- or z-axis. When the outflow is directed along the x axis, the accretion rate and acceleration are still irregular. The maximum outflow cases (f = 1) eject all gas accreted and shut off the accretion; for this reason they are not shown in the plots. Even with the maximum outflow efficiency, the acceleration decreases but still does not reverse direction.

Figure 13.

Varying the outflow efficiency f for simulations with outflow in the z-direction with u = 0.1 (V_w = 3 × 10³ km s⁻¹). Upper panel: the mass accretion rate. Lower panel: acceleration of the star in the x-direction. The dashed lines show the average value for t > 30t_B. The f = 1 line is not shown in the upper panel since all the accreted gas is ejected leaving a zero accretion rate.

Figure 14.

Varying the outflow efficiency f for simulations with outflow in the x-direction with u = 0.1 (V_w = 3 × 10³ km s⁻¹). Upper panel: the mass accretion rate. Lower panel: acceleration of the star in the x-direction. The dashed lines show the average value for t > 30t_B.The f = 1 line is not shown in the upper panel since all the accreted gas is ejected leaving a zero accretion rate.

4.3 Dynamical friction

In the previous section, we have shown that the existence of outflow will reduce total acceleration the compact object experiences. However, unlike the isotropic case where accretion is totally shut off by the outflow, the ambient gas keeps accreting to the central object. In our simulation where we set up a background with constant gas density, the overdense wake behind the object extends all the way to the boundary of the computational box. Due to the long-range nature of gravity, the value of total gravitational drag from the overdense wake increases with distance. The net gravitational acceleration contains two contributions: one from the overdense wake, which diverges with distance; and the other from the underdense region swept by the jet, which points to the opposite direction but converges with increasing box size.

To demonstrate this point, Fig. 15 shows the amount of gravitational acceleration a_G from gas inside a sphere with radius r from the object for different cases. When the outflow is absent, the analytical dynamical friction is proportional to the logarithm a_G ∝ ln Λ. As seen in Fig. 15, the blue line increases logarithmically with the radius r, as predicted by the analytical expression.

Figure 15.

Gravitational acceleration (not including accretion) in the x-direction computed from gas within a sphere with radius r. The acceleration is the average value for t > 0.01 yr. Upper panel: the jet is along the z-axis, perpendicular to the flow. Lower panel: the jet is along x-axis, parallel to the flow.

The same trend of increasing gravitational drag is observed for both cases where the jet is along the z- or x-axis (perpendicular or parallel to the flow, respectively). For large r, a_G grows logarithmically, signifying the domination of the gravitational drag from the overdense wake. The effect of outflow is more significant for smaller r, especially for powerful jet with small u (large V_w) and high efficiency f. We observe that the gravitational drag from a small sphere of gas r ≤ 4R_B points to the same direction as the motion of the object when the jet is along z-axis. For the cases where the jet is along x-axis, the gravitational acceleration are all positive. However, there is significant variability, as shown in Fig. 16. The variable acceleration does sometimes become negative of order 0.2a_DF, comparable with the value observed when the jet is along z axis. When the f is maximal, even though the average value of a_G does not change much, but a_G oscillates more wildly and gives rise to larger negative transient values.

Figure 16.

Evolution of gravitational acceleration in the x-direction computed from gas within a sphere of radius r = 5 × 10¹¹ cm for the case where the jet is along x- axis.

With these results in hand, we reconsider the isotropic wind case considered in Section 3. We note that, at sufficiently large distances, the positive force from the gravitational wake should overwhelm the negative dynamical friction observed the simulation presented in Section 3. This effect was not included in the analysis by Gruzinov et al. (2019), as this was based on Wilkin's (1996) mass distribution.

However, the underlying assumption of a uniform background also breaks down on large scales, because self-gravity would cause collapse on a time ∼(Gρ)^−1/2. (However this cannot occur in our simulations, thanks to the periodic nature of the self-gravity.) For the background fluid to be hydrostatic, therefore, its initial pressure scale height can be at most ∼c_s(Gρ)^−1/2. Significantly, in the astrophysical situations considered in Gruzinov et al. (2019), the spatial scale of the system was not that much greater than the Bondi radius. The finite size is clearly important both for the stellar mass black hole inspiraling through the common envelope, as well as for the stellar or intermediate-mass black hole embedded inside a thin active galactic nucleus (AGN) disc. The sign of the dynamical friction in these systems cannot be inferred from the numerical simulations in this paper, and requires a separate investigation.

5 DISCUSSION

We have performed hydrodynamical simulations of a compact object with an outflow moving through ambient gas of uniform density, focusing on the effects of the outflow and accretion on the dynamical friction. We considered the case of a constant, isotropic wind, the case of no outflow, and the case of a beamed jet with mass flux proportional to the accretion rate (either perpendicular or parallel to the flow direction). Our main results can be summarized as follows:

• Our simulations with strong isotropic winds lead to negative dynamical friction, as predicted by Gruzinov et al. (2019). The position of the bow shock and the strength of the negative dynamical friction are in agreement with analytical estimates. We anticipate that the net force is positive on sufficiently large scales, due to the gravitational wake; however the assumption of uniform density is inconsistent at large scales.

• Feedback from accretion-powered jets from the compact object is found to significantly reduce the accretion rate and the strength of dynamical friction in the BHL picture. The net force was never negative in our simulations when gravity of all gas from a large simulation domain was accounted for. However, the sign was reversed in some simulations, when gravity of the gas only within several Bondi radii was accounted for. This implies that the sign of the frictional force on compact objects with outflows needs to be treated with some caution, especially in dense environments.

We focus on the cases where the outflow is very powerful and the dynamical effect is most visible. As the outflow power is reduced, either with smaller velocity or smaller f, the dynamical effect will also be reduced and asymptote to the case with no outflow as seen in figures presented in Section 4.

Our simulations employed the simplest setup of a uniform background that extends to large distances (formally to infinity). Self-gravity of the uniform ambient gas was ignored. In physical systems, finite-size effects, density or velocity gradients, and self-gravity may all be important. Further complexities arise when the accreting object is itself orbiting within a disc of gas (in an AGN or protoplanetary disc), as in that situation forces are exerted at resonances whose position depends upon on the disc rotation profile. Our simulations should therefore be interpreted as local models (e.g. in the common envelope) that are approximately valid on scales where density gradients and other effects can be neglected. Though an outflow generally reduces dynamical friction in our simulations, whether dynamical friction can become negative depends on the competition between the negative effect of the outflow (primarily at small distances) and the overdense gravitational wake (at the largest distances). For a specific system, ad hoc simulations are required to determine the exact matter distribution at large distances.

In our simulations with accretion-powered jets, we have explored both cases where the jet is perpendicular and parallel to the object’s direction of motion. The perpendicular case shows a stable accretion and acceleration while the parallel case is more irregular and the accretion is strongly intermittent. In reality, the direction of the jet can be more irregular as the relative orientation between spin axis of the compact object and its velocity vector in the ambient gas’ frame can change during its motion. Irregularities in accretion rate and acceleration should be general features when the accretion-powered outflow is present.

In this paper, we assume the compact will accrete matter at the Bondi accretion rate |$\dot{M}_B$| if no outflow is active. A potential worry is that, in all our simulations, |$\dot{M}_B$| exceeds the Eddington luminosity, making the accretion hypercritical. We note, however, that the Eddington limit becomes irrelevant at sufficiently high accretion rate as the photons are trapped and advected inward with the flow (Rees 1978; Begelman 1979; Blondin 1986). As an example, Soria et al. (2014) report an accreting black hole in M83 with super-Eddington kinetic power. Bethe & Brown (1998) suggest that |$\dot{M}_B$| is the appropriate accretion rate for compacts in the common envelope phase. However, population synthesis of binary compact object systems (Belczynski, Kalogera & Bulik 2002; Belczynski et al. 2007) indicates that accretion at |$\dot{M}_B$| would cause neutron stars to collapse to black holes, which would be incompatible with observations. Hence, the true accretion rate must be reduced (Ricker & Taam 2008), and the exact process is still not well understood. Our simulations show that the dynamical effects of outflow can reduce the true accretion rate below |$\dot{M}_B$| which might serve as a partial solution to this problem.

In this work, we did not model the outflow launching mechanism in any detail and rather employed a phenomenological model. For a compact object of stellar mass, the outflow launching mechanism acts on scales much smaller than the dynamical scales that we considered, making simultaneous consideration of outflow launching and dynamical friction highly challenging. How the large- and small-scale dynamics are coupled and how the large-scale dynamics can change the disc formation and jet launching remain interesting questions to be explored.

ACKNOWLEDGEMENTS

We thank Andrei Gruzinov for useful discussions. PC is supported by the NASA ATP program through NASA grant NNH17ZDA001N-ATP and the Simons Foundation. We also use the yt software platform for the analysis of the data and generation of plots in this work (Turk et al. 2011). Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Colleges and Universities. The Flatiron Institute is supported by the Simons Foundation.

Footnotes

1

There are a variety of definitions of in the literature, but all lead to accretion rates that are essentially the same for cases in this paper.

2

In simulations, we take b_max to be the half box size and b_min to be the smoothing length of gravity.

3

We refer to a force opposite to the direction of motion as a positive dynamical friction force, and vice versa. Because the star experiences a wind from negative x, negative friction forces point along −x.

4

We do not study the effect of varying the accretor size in this paper. Though previous 2D simulations show that flow instabilities will develop for a sufficiently small accretor (Fryxell & Taam 1988; Blondin & Pope 2009), that effect is less significant in 3D simulations (Ruffert 1997, 1999; Blondin & Raymer 2012).

5

However, we note that equation (11) is based on a mass distribution derived by neglecting gravity; on sufficiently large scales we should expect the numerically derived force to eventually deviate from it – see Section 4.3.