Mechanical feedback effects on primordial black hole accretion

Dark matter may consist, at least partially, of primordial black holes formed during the radiation-dominated era. The radiation produced by accretion onto primordial black holes leaves characteristic signatures on the properties of the medium at high redshift, before and after Hydrogen recombination. Therefore, reliable modelling of accretion onto these objects is required to obtain robust constraints on their abundance. We investigate the effect of mechanical feedback, i.e. the impact of outflows (winds and/or jets) on the medium, on primordial black hole accretion, and thereby on the associated radiation. Using analytical and numerical calculations, we study for the first time whether outflows can reduce the accretion rate of primordial black holes with masses similar to those detected by the LIGO-Virgo collaboration. Despite the complexity of the accretion rate evolution, mechanical feedback is able to significantly reduce the primordial black hole accretion rate, at least by an order of magnitude, when outflows are aligned with the motion of the compact object. If the outflow is perpendicular to the direction of motion, the effect is less important but still non-negligible. Outflows from primordial black holes, even rather weak ones, can significantly decrease the accretion rate, effectively weakening abundance constraints on these objects. Our results motivate further numerical simulations with a more realistic setup, which would yield more precise quantitative predictions.

The characterization of accreting PBHs requires a number of assumptions that have not been fully tested yet. The robustness of PBH abundance constraints, coming from CMB, 21-cm signal, and galactic emission, depends significantly on those assumptions. In particular, one important aspect of PBH accretion-emission physics that has been little discussed in the literature concerns the indirect im-pact on accretion of outflows, winds and/or jets, produced by the PBHs themselves. This effect, known as mechanical feedback, consists of a reduction of gas accretion caused by outflows sweeping the medium away, as this leaves only diluted and hot gas to accrete. Mechanical feedback is in fact an important astrophysical process in general, and has been discussed for instance in the context of heavy stellarmass BHs as those detected by LIGO (Ioka et al. 2017), supermassive BHs growth (see e.g. Levinson & Nakar 2018;Zeilig-Hess et al. 2019, and references therein), and other astrophysical objects such as stars, galaxies and clusters (e.g. Soker 2016; Gruzinov et al. 2019;Li et al. 2019).
The aim of the present work is to estimate for the first time, using analytical and numerical tools, the potential reduction of the accretion rate, hence in the associated radiation, caused by the production of outflows in a moving PBH. Li et al. (2019) already tentatively explored through a numerical study the impact of mechanical feedback on Bondi-Hoyle-Lyttleton accretion in different astrophysical contexts. Here we consider the PBH case due to its far-reaching astrophysical and cosmological consequences. Thus, we focus on the case in which a PBH moves with a supersonic speed through an almost homogeneous medium, similar to that of the Universe at high redshift, when nonlinear structures had not formed yet, both before and after Hydrogen recombination. Therefore the results of this work can be applied both to existing CMB constraints and to future 21-cm constraints, but also have relevance for PBH Article number, page 1 of 10 arXiv:2004.11224v2 [astro-ph.CO] 24 Apr 2020 A&A proofs: manuscript no. vbosch (and isolated BHs in general) detectability in the presentday Universe.
The paper is structured as follows: first we provide an analytical characterization of the impact of outflows on the accretion process (Sections 2 and 3). Then, we present the results of numerical simulations to complement the analytical characterization (Sect. 4). Finally, we conclude with a summary and a discussion of our findings (Sect. 5). Throughout this work, we adopt the convention Q b = Q/10 b , where Q is any physical quantity measured in cgs units 1 , unless indicated otherwise. Also, the symbol ∼ indicates an order of magnitude estimate, whereas the symbol ≈ is used when we want to provide an approximate numerical estimate of a given quantity.

Accretion
We consider a PBH of mass M PBH that moves with supersonic speed v PBH through an almost homogeneous medium with density ρ m = m n m , number density n m , and average particle mass m ∼ m H ≈ 1.7 × 10 −24 g, that of the Hydrogen atom mass. We assume as reference case a PBH with mass M PBH = 30 M moving with a speed of v PBH = 3×10 6 cm s −1 in a medium with number density n m = 1 cm −3 (i.e. M PBH,1.5 ≈ 1, v PBH,6.5 ≈ 1, n m,0 = 1).
As proposed in Ali-Haïmoud & Kamionkowski (2017), if PBHs are the dark matter, their typical speed is approximately constant before Hydrogen recombination, with v PBH ≈ 3 × 10 6 cm s −1 (Dvorkin et al. 2014), and v PBH ∝ (1+z) afterwards (Tseliakhovich & Hirata 2010), where z is the redshift. Since PBHs move supersonically with respect to baryons until relatively low redshift (z ≈ 20), our set-up roughly captures both pre-and (early) post-recombination scenarios. The reference value taken for the number density of the medium, n m = 1 cm −3 , corresponds to the baryon number density n b (z) = n b0 (1 + z) 3 at redshift z ≈ 150, where n b0 ∼ Ω b0 ρ 0c /m H ≈ 3 × 10 −7 cm −3 is the presentday baryon number density, Ω b0 ≈ 0.05 the baryon fractional abundance at present time, and ρ 0c ≈ 10 −29 g cm −3 the present-day critical density. 2 A supersonic PBH accretes approximately at the Bondi-Hoyle-Lyttleton accretion rate (Hoyle & Lyttleton 1939;Bondi & Hoyle 1944): where the effective accretion radius r acc , which defines the PBH sphere of influence radius, is given by: (2) 1 Useful conversion factors between different systems of units: 1 g s −1 ≈ 10 −26 M yr −1 ≈ 10 21 erg s −1 /c 2 , where the speed of light c in cgs units is c 2 ≈ 10 21 cm 2 s −2 ; the gravitational constant in cgs units is G ≈ 6.7 × 10 −8 cm 3 g −1 s −2 ; 1 cm ≈ 3 × 10 −19 pc; and 1 s ≈ 3 × 10 −8 yr. 2 We note that, in the adiabatic case of mechanical feedback, our conclusions will not depend on nm, and thus its actual value is not of much relevance here. Fig. 1: Sketch of the no-outflow case, where the PBH is moving supersonically through a medium, and it is accreting material from its sphere of influence through a gas column forming opposite to the direction of motion. Sketch not to scale.
In case the PBH moves with subsonic speed, one has to substitute v PBH in Eq. (2) with the sound speed of the medium (Bondi 1952). In that situation, accretion would be closer to the spherical case, at least on scales r acc (but not necessarily in the close vicinity of the PBH; Waters et al. 2019). This would be the case for instance at redshift z 20, but we do not investigate that epoch here due to its high complexity, as it coincides with the formation of the first structures.
Processes in the PBH vicinity may reduce the accretion rate in that region. Typically, they are accounted for by defining the real accretion rate asṀ = λṀ PBH , wherė M PBH is the accretion rate defined in Eq. (1), and the parameter λ 1 encloses the physical effects of pressure, viscosity, radiation feedback and other processes on scales close to the PBH gravitational radius (r g = GM PBH /c 2 ≈ 5 × 10 6 cm for M PBH = 30 M ). In fact, λ captures the presence of gravitationally unbounded flows of any type. In particular, the production of fast winds or jets can significantly reduce the accretion rate on those scales (i.e. λ 1; see e.g. Tsuna et al. 2018;Bu & Yang 2019, in the context of isolated stellar-mass and supermassive BHs). Here, since we are interested only on mechanical feedback on scales comparable to r acc , any process that reverts accretion close to the PBH is phenomenologically treated as outflow on larger scales. Thus, our reference value for the accretion rate isṀ PBH , i.e. we take λ = 1.
Relatively far from the PBH, accretion takes place through a column that forms behind the accretor, as schematically drawn in Fig. 1. In the case of a perfectly homogeneous medium, accretion is axisymmetric, and matter is channeled by this gravitationally bounded gas column straight to the accretor. However, it is expected that small perturbations, certainly present in a cosmological context, will lead to a break of the axisymmetry of the accretion process, which can be then enhanced by instability growth.
This can lead to the accretion of some angular momentum (see e.g. Foglizzo et al. 2005;Lora-Clavijo et al. 2015; see also Waters et al. 2019 for the spherical case). In fact, density inhomogeneities accreted with different impact parameters will already provide some angular momentum. In these contexts, even if the average angular momentum accreted in the long-term is very small, it is possible to form a short-lived small disc close to the PBH (see e.g. Agol & Kamionkowski 2002, and references therein).
To further characterize accretion in the vicinity of the PBH, it is informative to compare the BH accretion luminosity to the Eddington luminosity. Accounting only for Thomson scattering on ionized Hydrogen, the Eddington luminosity is given by L Edd ≈ 3.8 × 10 39 M PBH,1.5 erg s −1 . Defining here the Eddington accretion rate asṀ Edd = 10L Edd /c 2 , we obtain that the ratio of the PBH accretion rate to the Eddington one is given bẏ which shows that, for a very broad range of parameter values,Ṁ PBH /Ṁ Edd 0.1. Such a low Eddington accretion fraction implies that in our scenario the accreted gas will not cool efficiently on scales r g , leading to the formation of a faint and hot (virialized) accretion structure close to the PBH (see e.g. Shakura & Sunyaev 1973). As we will see, such a configuration is prone to outflow formation.
Observational cosmology is unable to establish whether inhomogeneities lead to net angular momentum accretion because of the small spatial scales involved. The key quantities that regulate angular momentum accretion are amplitude and spatial scale of the inhomogeneities. Regarding the former, we note that the development of instabilities in the inner accretion structure can help to amplify rather smooth inhomogeneities. Concerning the latter, by characterizing the inhomogeneity scale as l inh , one can estimate the duration of net angular momentum accretion: t l ∼ l inh /v PBH . For inhomogeneities to be important, it is required that: (i) t l r g /c, i.e. t l is larger than the wind formation time scale; and (ii) t l r acc /v PBH , i.e. t l is smaller than the PBH sphere of influence crossing time. Otherwise, (i) the wind will not have time to form, or (ii) the gradients may be too small to accrete a non-negligible amount of net angular momentum. Given the random nature of the inhomogeneities, both in amplitude and in l ihn , it is expected that the outflow will change orientation on a time ∼ t l . Given that the evolution timescale of the whole system is expected to be longer than any relevant t l (see below), the outflow could effectively act, on average, as a very broad wind at large scales.

Outflows
A bipolar outflow (wind and/or jet) from an accreting PBH can affect its own energetics, and the radiation from accretion, by sweeping away medium gas on scales ∼ r acc . Nevertheless, the launch of such an outflow requires: (i) a strong magnetic field of specific geometry; and/or (ii) an excess of thermal energy in some of the accreted gas, i.e. positive energy gas that turns into a wind.
Regarding (i), one must note first that the magnetic field in the accreted primordial gas is expected to be very small, although the Biermann battery mechanism may alleviate this problem in the accretion flow, and the magnetorotational instability could produce an enhancement of magnetization and turbulence, thereby generating the viscosity needed for accretion (see e.g. Safarzadeh 2018) 3 . However, the magnetic field seed generated by the Biermann battery is toroidal, whereas jet launching requires some vertical magnetic flux (see e.g. the discussion in Sotomayor Checa & Romero 2019). Nevertheless, if enough vertical flux is present close to the PBH, the inner accretion region can efficiently launch a collimated, supersonic bipolar outflow (Blandford & Payne 1982). In addition, if the PBH is spinning fast, relativistic, highly collimated jets can be also produced in the PBH ergosphere through the Blandford-Znajek (Blandford & Znajek 1977) mechanism (even if accretion is quasi-spherical, as discussed e.g. in ).
In the case of (ii), regardless of the presence of jets, a weakly accreting, geometrically thick disc can form winds through thermal pressure gradients (see e.g. Sadowski et al. 2016, where jets are also shown to form if the BH is spinning). Such (less collimated) outflows, produced by an excess of thermal energy in the accretion structure, are less dependent on the details of the magnetic field. Thus, to be conservative, hereafter we assume that the outflow has a thermal origin, and carries only a small fraction of the rate of accreted energy ∼Ṁ PBH c 2 . The power of each constituent of the bipolar supersonic outflow can be characterized as: Alternatively, one can write in the non-relativistic and the highly relativistic (cold outflow) case, respectively, where v w is the outflow speed,Ṁ ejec the total ejected mass, and γ w the outflow Lorentz factor. Equations (4) and (5) implẏ which shows that v w 6 × 10 8 1/2 −4 cm s −1 to ensure thaṫ M ejec <Ṁ PBH (i.e. large values require highly relativistic outflows).
One can also estimate a minimum for the outflow to escape the vicinity of the PBH: The outflow propagation should proceed roughly unimpeded if accretion is rather equatorial, away from the outflow close to the PBH. However, beyond some distance, hard to ascertain but prescribed here by the parametric length h sph , the inflowing gas could be quasi-isotropic. If such a region were present within the PBH sphere of influence, the supersonic outflow should have enough ram pressure to overcome that of the inflowing gas, if distances > r acc are to be reached. One can estimate when it is the case assuming that accretion is roughly spherical above h sph from the PBH, that the gas flows inwards at free fall velocity, and then deriving a constraint on by equating the outflow and medium ram pressures, p w = p m at h sph : (9) for the non-relativistic and the highly relativistic case, respectively; χ is the outflow half opening angle. Thus, the outflow will escape if or for the non-relativistic and the highly relativistic case, respectively, where tan χ has been approximated as χ (valid for χ 1).

Mechanical feedback: analytical estimates
In this section, we characterize analytically the interaction of the outflow with the surrounding medium in its early (Sect. 3.1) and long-term (Sects. 3.2 and 3.3) stages of evolution. We discuss the effect of the relative orientation between the outflow and the PBH velocity, corresponding to different values of the angle θ shown in Fig. 2. We investigate the cases when the outflow and the PBH velocities are orthogonal (θ = π/2; Sect. 3.2), and when they are aligned (θ = 0; Sect. 3.3). In the long-term evolution, the value of the half-opening angle χ is expected to be also important for mechanical feedback; we discuss two qualitatively distinct cases: collimated outflows with χ 1 (see Sects. 3.2 and 3.3) versus non-collimated outflows with χ ∼ 1 (see Sect. 3.3). In this section we consider as reference case a conical outflow with half-opening angle χ = 10 −1 , velocity of the outflow v w = 10 9 cm s −1 , and luminosity L w = 10 30 erg s −1 , corresponding to ≈ 10 −4 (see Eq. 4) for our choice of parameters. Nevertheless note that our predictions will not be very sensitive to the actual parameter values (with the clear exception of χ) as long as L w is large enough to allow the outflow front to escape the PBH sphere of influence.

Early interaction
After leaving the region close to the compact object, the front of the supersonic outflow is eventually slowed down by the ram pressure of the gas surrounding the PBH. Since the medium sound speed is significantly smaller than v w , the medium ram pressure in the outflow reference frame is The ram pressure of the outflow front in the medium frame is: in the non-relativistic case, where S = πh 2 tan 2 χ is the outflow front surface at the distance h from the PBH. For typical values of the parameters we have that the outflow and medium pressures become comparable at a height: where tan χ ∼ χ. If the outflow velocity were highly relativistic, then p m ≈ ρ m γ 2 w v 2 w and p w ≈ L w /(Sv w ), turning h into h rel = h/(γ w √ 2). Approximately at the distance h from the PBH, where pressure equilibrium is reached, the front of the supersonic outflow significantly slows down and the material reaching the front gets shocked. Then, this shocked material is directed sideways and then backwards, and inflates a bubble of hot adiabatic gas surrounded by a shell of shocked medium. The outflow is initially freely expanding, but at some point the bubble pressure balances the lateral pressure of the outflow, point at which the latter acquires a cylindrical geometry. This leads to the self-similar evolution of the bubble, in which the bubble dynamics depends on its energy content and the accumulated medium mass (see e.g. Alexander 2006, in an extragalactic context). The bubble evolution is thus not very sensitive to how collimated or relativistic the initial outflow was. The described process is illustrated in Fig. 2.
The self-similar bubble expansion can be approximated as adiabatic and symmetric, fed by a steady source of energy; thus, its evolution can be approximately characterized by the Sedov-Taylor equation in the continuous regime (see e.g. Blandford & McKee 1976 and references therein; see also Kaiser &Alexander 1997 andAlexander 2006 in the context of extragalactic jets), according to which the bubble radius R b evolves in an approximate self-similar fashion Article number, page 4 of 10 V. Bosch-Ramon and N. Bellomo: Outflows and accretion in primordial black holes as: with forward velocity: Thus, for some time the bubble is expanding forward significantly faster than the PBH moves, for v PBH = 3 × 10 6 cm s −1 . Here, for simplicity, we assume that the bubble is approximately spherical, although the lateral expansion can be a factor of a few slower (e.g. Kaiser & Alexander 1997), which should be accounted for in a more detailed treatment of the problem. If the pressure equilibrium distance is larger than the effective accretion radius, i.e. h > r acc , the bubble will start to form after a time t h ∼ h/v w , and once the outflow front has already left the PBH sphere of influence. In this case, the time required by the front to cross the sphere of influence t cr roughly corresponds to the time needed to reach a distance r acc from the PBH: t cr ∼ 10 6 r acc,15 /v w,9 s (< t h ). Viceversa, if h < r acc , the crossing time will be t cr ∼ r acc /v b (t cr ) ≈ 10 7 r 5/3 acc,15 n 1/3 m,0 L −1/3 w,30 s r acc /v w . In either case, after a time max[t h , t cr ], the bubble filled with hot gas will engulf the PBH sphere of influence and accretion is expected to be affected. Only if v b becomes v PBH before the bubble reaches out of the PBH sphere of influence (i.e. an extreme case plus the condition h < r acc ), bubble growth will be overcome by gravity, and the outflow will not be able to prevent accretion.

Quasi-perpendicular outflow
Once outside the PBH sphere of influence, the bubble slows down until its expansion velocity and the PBH velocity become comparable, i.e. v b ∼ v PBH . At that moment and for θ ∼ 1, the bubble is already far from symmetric because the medium ram pressure is mostly acting on one side of the bubble. The typical timescale of the previous axisymmetric expansion t ax can be derived from Eq. (15): and the bubble size at time t ax is: ≈ 5 × 10 16 v PBH,6.5 t ax,10 cm.
Therefore, for the system under consideration, we have that after t ∼ t ax ≈ 400 yr the bubble will reach a typical size of ∼ 60 r acc (see Eqs. 2 and 17). At that point, the bubble is already getting compressed in the direction of the PBH motion, and tends to expand in the opposite direction, as shown in Fig. 3 (see also, e.g., Yoon et al. 2011 for simulations of jets from a BH microquasar moving in the The strong bubble asymmetry makes the two constituents of the bipolar outflow deflect away from the PBH motion direction. The typical scale d at which the outflow deflection is significant can be derived from equating the medium ram pressure in the PBH reference frame, p m ≈ ρ m v 2 PBH , with the outflow ram pressure, p w ≈ 2L w /v w πd 2 χ 2 in the non-relativistic case: where we assume that χ 1. In the highly relativistic case, for which p w ≈ L w /cπ(d rel ) 2 χ 2 , this distance turns into d rel = d/ √ 2. The deflection distance d turns out to be much larger than the PBH sphere of influence; for the parameters considered in this work: d ∼ 70 r acc . This scale also characterizes the typical size of the section of the interaction structure perpendicular to the PBH motion.
For θ = π/2 (and χ 1), free sidewards expansion of the outflow propagating within the asymmetric bubble is stopped on the side on which the medium ram pressure acts, on a distance d r < d from the PBH. Assuming for simplicity that the outflow is very cold, outflow asymmetric reconfinement occurs at the distance where the medium and outflow lateral ram pressures are equal, ρ m v 2 PBH = 2L w /πd 2 r v w in the non-relativistic case: In the highly relativistic case, where ρ m v 2 PBH = L w /γ 2 w cπ(d rel r ) 2 , this distance becomes d rel r = d r /(γ w √ 2). Within d r , the outflow is expected to be ballistic, at ∼ d r it starts to deflect, and at d, it has significantly deviated. The tail described above is in fact made of shocked deflected outflow heavily mixed with shocked medium, as the outflow-medium contact discontinuity is prone to hydrodynamical instability growth (see e.g. Sect. 4).
As shown by Li et al. (2019), for θ = π/2 the surrounding medium gas can penetrate into the sphere of influence and reach the PBH on the directions perpendicular to the outflow, because the bubble pressure is the lowest there. Thus, accretion may only be moderately affected by mechanical feedback. Nevertheless we note that for χ 1 and θ < π/2 (but still ∼ 1), the outflow can suffer a stronger, oblique shock that can widen the outflow front. This situation might turn into a case where we effectively have χ ∼ θ, which is closer to the quasi-parallel case reported in Sect. 3.3.
So far we have assumed that the evolution of the outflow-medium interaction structure was adiabatic. However, according to the cooling rates for low metallicities in Sutherland & Dopita (1993) and the timescale of the problem ( t ax ) close the recombination epoch the medium was dense enough (n m (z rec ≈ 1100) ≈ 400 cm −3 ) to make the bubble-driven shock in the medium radiative. This affects the dynamics of the bubble; in particular, it enhances the growth of instabilities at the outflow-medium contact discontinuity. A fully disrupted contact discontinuity would make an analytical treatment of the studied scenario more speculative. At this stage, we leave for future work a more detailed analysis of the radiative bubble case.

Quasi-parallel outflow
When the outflow is roughly aligned with the direction of motion, accretion onto the PBH sphere of influence is expected to be almost halted once the hot and diluted bubble gas fully occupies that region, keeping out the external gas (as shown in Sect. 4, and also in Li et al. 2019). As already mentioned, this could occur also for non-negligible values of θ if the outflow is non-collimated, i.e. χ ∼ 1. Under these circumstances, significant outflow production can last only until the remaining material within the sphere of influence has been accreted.
The typical timescale of accretion of the remaining medium gas is of order of the free-fall timescale at r acc : After a time t cr + t acc from the formation of the outflow, only diluted hot matter should be present around the PBH and accretion should be very weak. The bubble size and expansion speed at this stage are given by Eqs. (14) and (15), and they read as and v b (t cr + t acc ) ≈ 2 × 10 7 R b,16 t −1 9 cm s −1 , respectively. Therefore, even after having expanded up to ∼ 20 r acc , the bubble expansion speed is still significantly faster than the PBH, for v PBH = 3 × 10 6 cm s −1 . When both accretion and ejection stop, energy injection into the bubble stops as well. At that point, the outflow has injected an energy of Once the energy injection stops, for t (t cr + t acc ) the bubble will expand adiabatically with size and velocity approximated as (see e.g. Blandford & McKee 1976, for the adiabatic and impulsive regime) and v b (t) ∼ 2 5 respectively. The time needed for the bubble velocity v b to reach the PBH velocity v PBH reads as: and the bubble final size after a time t cr + t acc + t end ∼ t end is: ,9.5 cm, (27) which is somewhat smaller than the value obtained in Eq. (18) for the quasi-perpendicular case, but still much larger than the PBH sphere of influence.
Once v b gets v PBH , the ram pressure of the medium becomes dominant over that of the shocked outflow, and the material starts re-approaching the PBH as the latter moves, reaching it in an accretion recovery time As (t rec + t end ) (t cr + t acc ), the stage in which the PBH neither accretes nor ejects largely dominates the duration of the whole on-off accretion cycle, and so the time-average accretion rate should be much lower than the Bondi-Hoyle-Lyttleton one. For the typical values derived here, we find an average reduction inṀ PBH by a factor (t cr +t acc )/(t rec + t end ) ∼ 0.1 (see Sect. 4), but lower and higher values are also possible given the large plausible parameter space.
The maximum, temporary section radius of the outflowmedium interaction structure in the quasi-parallel case can be identified with the bubble radius R b (t end ). If otherwise accretion and outflow production were constant, the structure section radius would be also constant and characterized by Eq. (18).

Mechanical feedback: numerical calculations
In Sect. 3.2 we have seen that a perpendicular, well collimated outflow may leave the PBH sphere of influence without strongly affecting the gas accretion from directions away from the outflow, allowing the gas to reach the PBH. On the other hand, even in the perpendicular case, accretion could be reduced if high-pressure, shocked outflow material were more homogeneously distributed inside the PBH sphere of influence, more effectively preventing the external medium from entering that region. As noted, this could be the case for more oblique and/or broad outflows, but even collimated (χ 1) perpendicular outflows may in principle divert enough momentum sideways due to instability growth closer to the PBH, leading to a situation closer to that discussed in Sect. 3.3. Unfortunately, finding the value range of θ and χ for which accretion is inhibited by mechanical feedback requires a number of expensive 3-dimensional (3D) hydrodynamical simulations, as they should include spatial scales both r acc and d; for typical parameter values this spans several orders of magnitude. Adaptive mesh refinement can be useful to reduce the computational time, but note that highly turbulent flows, as expected in the studied scenario, are sometimes challenging for this technique (see e.g. Bosch-Ramon et al. 2015).
The computing time needed for a simulation is proportional to the number of computing cells of the grid, times the number of simulation time steps: T sim ∝ N cells × N steps , where N cells ∝ (d/r acc ) D ∝ (v w /v PBH ) D/2 , with D being the dimensionality, and N steps ∝ t phys /∆t phys , where t phys ∝ d/v PBH is the total simulated time, and ∆t sim ∝ d/v w the simulated time step. This yields T sim ∝ (v w /v PBH ) D/2+1 . To make the problem more tractable, one can choose parameter values to reduce T sim , although one still requires that d/r acc 1 to capture at least qualitatively the main features of the system. This additional constraint makes a 3D simulation quite demanding. On the other hand, a 2D simulation, which is realistic in the parallel outflow case because of axisymmetry, becomes more affordable.
For this work, as a test of the case in which accretion is strongly reduced due to mechanical feedback, we carried out axisymmetric hydrodynamical simulations for θ = 0. We adopted values for v PBH and v w that allow us to probe a moderately large range of spatial scales d/r acc , and that reduce significantly N steps . An intermediate value of the outflow half-opening angle was obtained by adopting a Mach number of M ch = 2.5, which yields shortly after injection, already in the supersonic regime, a χ ≈ 1/M ch = 0.4; but we note that trials with χ = 0.1 (not shown here) yielded similar results.

Properties of the simulation
We performed axisymmetric hydrodynamical simulations of the interaction between a collimated outflow with v w = 10 9 cm s −1 and initial M ch = 2.5, and a moving medium with v PBH = 5 × 10 7 cm s −1 and Mach number of 3. The medium density was fixed to ρ m = 1.7 × 10 −24 cm −3 (i.e. n m ≈ 1 cm −3 ), and the power of each constituent of the bipolar outflow to L w = 0.025Ṁ PBH v 2 w , corresponding to ≈ 10 −4 , which is the reference value used for our analytical estimates above. The grid size was chosen such that the interaction structure remained in the computational grid for the duration of the simulation. We note that exploratory trials done with a higher -value suggested that results should not change significantly.
The cell size of the simulation was fixed to 10 11 cm. The size of the grid in the radialr-and the verticalẑ-directions are 250 cells (2.5 × 10 13 cm) and 500 cells (5 × 10 13 cm), respectively. The source of the gravitational potential is assumed to be point-like, with mass M PBH = 30 M , but the grid region where matter is accreted and leaves the grid (i.e. the accretor in the context of the numerical simulation) is modelled as a spherical region of low density and pressure, with a radius of 10 cells and centred at the origin. The outflow is injected in the ±ẑ-directions through a cylindrical inlet with a height and radius of 20 and 5 cells, respectively, and centred also at the origin, occupying thus a fraction of the accretion sphere. Any material crossing the accretor spherical surface but at the outflow inlet boundaries immediately disappears from the grid. Outside the outflow inlet and the accretor region, reflection conditions are imposed at the vertical axis (r = 0: left grid boundary), outflow conditions 4 at the upper and right grid boundaries, and inflow conditions (i.e. the medium) at the bottom grid boundary. The motion of the medium in the PBH reference frame is in the +ẑ-direction. The accretion rate is computed on the cells right outside the boundary of the accretor (i.e. excluding the outflow inlet). In this specific set-up the sphere of influence has approximately a radius of 30 cells (3 × 10 12 cm). Trials carried out with better resolution by a factor ≈ 1.9, and a slightly smaller accretor radius, ≈ 26% instead of 30%, yielded very similar results to those presented here.
The code used to solve the hydrodynamics equations is the same as the one used in de la Cita et al. (2017): third order in space (Mignone et al. 2005); second order in time; and using the Marquina flux formula (in the nonrelativistic case, Donat & Marquina 1996). Both the outflow and medium fluids were simplified as non-relativistic mono-atomic gases with adiabatic index 5/3. We neglect the effects of ionization and gas abundances, which are not relevant on the simulated scales ( r g ), on which gravity and ram pressure dominate the dynamics of an adiabatic flow. Any magnetic field should be affected by the flow evolution, but here we are assuming that the role of the magnetic field is purely subsidiary to the collisionless gas hydrodynamical approximation, and so of no dynamical importance. The accretor gravitational force was appropriately included in the hydrodynamical equations as a source term (see e.g. Toro 2009).

Results
We show in Fig. 5 the accretion rate of two cases: when the mechanical feedback mechanism due to the outflow is present and when it is not. When there is no feedback, the long-term accretion rate converges to a value within a factor of 2 fromṀ PBH ≈ 3 × 10 −17 M yr −1 predicted by Eq.
(1). On the other hand, when outflows are launched, our simulations indicate that matter reaches the PBH only indirectly, through filamentary structures related to hydrodynamical instabilities in the contact discontinuity between the outflow and the medium. Thus, accretion is on average one order of magnitude lower (Ṁ PBH ≈ 3×10 −18 M yr −1 ) than in the case without outflow, but presents strong fluctuations with respect to that average value. It is worth indicating the presence of filamentary material entering the grid from the top for t 15 t acc . This backflow is a numerical artifact associated to the boundary conditions and the presence of subsonic material in the upper boundary. The subsonic material in the boundary can induce the flow to bounce back, sending waves into the grid. However, we have checked with additional simulation trials that imposing a positive velocity at the boundary (i.e. the flow is forced to leave the grid), of modulus 5 × 10 7 cm s −1 , quantitatively yields very similar results. The simulation was run until it presented repetitive patterns of the flow, a sort of quasisteady state. Our results are quite similar to those obtained by Li et al. (2019) for θ = 0. We also present in Fig. 6 the density maps showing the evolution of the outflow structure for approximately 50 times the accretion dynamical timescale t acc (t acc ≈ 6 × 10 4 s, according to Eq. 20). The time of each snapshot in Fig. 6 is marked with a cross at its corresponding location in the accretion rate curve of Fig. 5. These density maps can be compared to the no-outflow simulation, run for a very similar amount of time, and shown in the bottom right panel of Fig. 6. Both Figs. 5 and 6 are illustrative of the strong non-linear nature of the outflow-medium interaction structure, and reminiscent of the intermittent nature of accretion predicted in Sect. 3.

A toy-model for θ = π/2
Despite the difficulty of studying the quasi-perpendicular case numerically, we still made a toy simulation in 2D of the case θ ∼ 1 by introducing an additional gravitational potential with a velocity dispersion σ, and with the medium initially at rest (as under the presence of a dark matter halo; see e.g. Zeilig-Hess et al. 2019, for a calculation in the supermassive black-hole case). We ran test simulations of that scenario in the context of a PBH of 30 M and σ = 30 km s −1 and of a BH of 10 8 M and σ = 300 km s −1 , and obtained similar results in both cases to those of Zeilig-Hess et al. (2019) for the second case, i.e. the accretion rate tended to a value a factor of a few smaller than the Bondi value without additional gravitational potential 5 . This result cannot be overstressed as it is not 3D, but shows the robustness of our numerical scheme. Moreover, the derived reduction in the accretion rate is similar to that found by Li et al. (2019) when computing a case in actual 3D similar to the one discussed in Sect. 3.2. This similarity is likely related to the fact that in both cases, 3D supersonic accretion with outflow and θ = π/2, and 2D subsonic accretion with outflow and adding a dark-matter-like potential, there is a characteristic medium velocity (v pbh and σ, respectively) leading to reduced but steady accretion on the directions away from the outflow.

Discussion and conclusions
The effect of mechanical feedback is expected to reduce accretion onto a moving massive object, but quantifying this effect is difficult. In this work, we show that the angle between the outflow direction and the PBH motion, θ, is likely to determine how far below the Bondi rate the PBH accretion is. A numerical study of a similar but more general scenario by Li et al. (2019) reaches results in line with those presented here in the context of PBH. The half-opening angle χ is also expected to play an important role, as broader outflows are expected to have larger solid angles for which accretion is halted. As our analytical estimates and simulations show, an outflow do not need to be strong in order to reduce accretion; it is enough that the outflow effectively heats and dilutes all the gas in the PBH surroundings on scales r acc .
Assessing the level of accretion as a function of θ and χ is challenging due to the complexity of the outflow-medium interaction, and its intrinsic 3D nature in general. Nevertheless, we can already predict that accretion is expected to be significantly smaller than Bondi-Hoyle-Lyttleton accretion, with the presence of strong fluctuations, at least for values of θ χ. For θ χ, analytical and numerical calculations seem to indicate that the reduction may still be non-negligible, but with accretion proceeding more steadily. Nevertheless, the larger the ratio v w /v PBH , the larger the region affected by the outflow, in which case there is more room for some redirection of outflow momentum and energy, caused for instance by instability growth on the outflow walls. This might lead to the formation of a more effective barrier for accretion on scales similar to r acc , even for θ ≈ π/2 and χ 1. In addition, θ may change with time, as explained at the end of Sect. 2.1, which could effectively lead to a time-averaged χ ∼ 1. For all this, knowing the effect of the time-averaged values of θ and χ of a PBH, and their statistical distribution in the PBH population(s), is crucial to estimate the effects of their accretion and related emission.
It is clear that the inclusion of mechanical feedback on a PBH population has a direct impact on PBH abundance constraints, as they depend on the energy injected by the PBH population into the primordial medium. This can be exemplified by the following. In the case of CMB constraints (Ricotti et al. 2008;Ali-Haïmoud & Kamionkowski 2017;Poulin et al. 2017;Bernal et al. 2017;Nakama et al. 2018), the injected energy density rate scales asė PBH ∝ f PBHṀ 2 PBH , where f PBH is the fraction of dark matter made up of PBHs. Fixingė PBH and assuming for simplicity the same M PBH andṀ PBH for all PBH, a decrease of an order of magnitude inṀ PBH would allow an increase of two orders of magnitude in f PBH . This would largely relax the present constraints on f PBH in the 10−100 M mass range. In addition to CMB constraints, constraints derived from galactic emission (Gaggero et al. 2017;Manshanden et al. 2019) could be also affected by mechanical feedback.
We have not discussed here the effects on the medium of the non-thermal radiation directly generated by the PBH outflows, because in our case outflows are rather weak, i.e.
1. Notwithstanding, outflows from PBHs may have ≈ 1 and high radiative efficiency, in which case their highenergy radiation could contribute to heat, excite, and ionize the medium in addition to direct accretion emission. This would make the constraints on PBH abundance tighter, so such an effect should be considered in the future. In fact, radio (Gaggero et al. 2017;Manshanden et al. 2019) and high-energy emission from outflows produced by PBHs, located in dense regions of our galaxy, may be important enough to be detectable at galactic distances. This idea has been already explored in  in the case of isolated stellar BHs. Fig. 6: Colour density maps of the outflow-medium interaction structure after (from left to right, and top to bottom) t ≈ 0, 5, 10; 15, 20, 25; 30, 35, 40; 45, 50 t acc . The bottom right density map corresponds to the case without outflow after t ≈ 50 t acc . The arrows illustrate the accreted gas trajectories. The injector-accretor region is located at (0, 0). The axis scales are normalized to r acc , and the colour scale is normalized to the medium density. The medium motion in the PBH reference frame is upwards.