{'fraction_non_alphanumeric': 0.03120971777801368, 'fraction_numerical': 0.01424000242641149, 'mean_word_length': 4.766681241801487, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The effect of outflows on CMB bounds from Primordial Black Hole accretion

Should Primordial Black Holes (PBHs) exist in nature, they would inevitably accrete baryonic matter in their vicinity. In turn, the consequent emission of high-energy radiation could affect the thermal history of the universe to an extent that can be probed with a number of cosmological observables such as the Cosmic Microwave Background (CMB) anisotropies. However, our understanding of the accretion and radiation emission processes in the context of PBHs is still in its infancy, and very large theoretical uncertainties affect the resulting constraints on the PBH abundance. Building on state-of-the-art literature, in this work we take a step towards the development of a more realistic picture of PBH accretion by accounting for the contribution of outflows. Specifically, we derive CMB-driven constraints on the PBH abundance for various accretion geometries, ionization models and mass distributions in absence and in presence of mechanical feedback and non-thermal emissions due to the outflows. As a result, we show that the presence of such outflows introduces an additional layer of uncertainty that needs to be taken into account when quoting cosmological constraints on the PBH abundance, with important consequences in particular in the LIGO-Virgo-KAGRA observational window.


Contents
1 Introduction The ΛCDM model has become the standard cosmological model due to its remarkably accurate predictions of a variety of observables, ranging from the Cosmic Microwave Background (CMB) radiation [1] to the formation and evolution of large scale structures [2][3][4]. 1 However, this phenomenological model does not provide any fundamental explanation for the nature of its "dark" components, i.e., cold Dark Matter (DM) and dark energy, despite the fact that they make up for the vast majority of the energy density of the universe today. Therefore, in order to find an origin for these fundamental, yet elusive, components of the universe, arguments invoking non-standard (beyond-ΛCDM) physics are often put forward. In particular, in the context of DM, Primordial Black Holes (PBHs) have regained interest in recent years as potential DM candidates [9][10][11] after the first LIGO detection of a Gravitational Wave (GW) signal emitted from a binary BH merger [12]. PBHs are BHs that formed in the very early universe, much before the appearance of the first stars. Their existence was first suggested in the late '60s [13][14][15][16] and they have since become one of the most popular candidates to make up for at least a fraction of the total DM content of the universe (see, e.g., refs. [17][18][19][20] for recent reviews). Moreover it has been suggested that PBHs could also be the seeds for the formation of intermediate mass [21] as well as super massive [22,23] BHs, which are believed to reside at the center of most galaxies. Other observational conundra that are possibly solved by the potential presence of PBHs are discussed in refs. [19,24].
Despite these arguments in support of the presence of PBHs, many of their properties, such as their formation mechanism, abundance or mass distribution, are still unknown [19,24,25]. As a consequence, existing constraints on the PBH abundance heavily rely on assumptions about the BH phenomenology. For instance, it is well known that PBHs, if existing, would inevitably accrete matter and convert a fraction of the energy released by the accretion into radiation. The resulting injection of high-energy photons into the cosmological photon bath would then affect the thermal history of the universe and, thus, observables that depend upon it, such as the CMB anisotropies. However, the impact of this energy injection on the CMB varies depending on the details of the accretion and energy injection mechanisms considered, for instance on the assumed geometry of the accretion [26][27][28], or on the properties of the environment surrounding the PBH [29,30]. Furthermore, the potential presence of magnetic fields and/or the excess of thermal energy in part of the accreted matter might lead to the formation of outflows, like winds or jets. Such outflows are expected to dilute the medium surrounding the BH reducing the total amount of accreted material, effect known as mechanical feedback (MF), and to be potential accelerators of non-thermal particles capable to effectively boost the total luminosity of the BH [31][32][33][34][35]. Overall, the resulting CMB constraints on the PBH abundance cover approximately the O(1 − 10 4 ) M mass range, with vastly different outcomes depending on the accretion and emission mechanisms [27,28,30,36]. These CMB anisotropy constraints are also complemented by bounds derived from CMB spectral distortions [27], galactic emission [37][38][39], interstellar gas heating [40], and the 21-cm line [41][42][43].
Yet, an accurate modelling of the accretion mechanism is fundamental to test the consistency between cosmological constraints and astrophysical hints in support of the PBH hypothesis. In fact, the range of PBH masses that would explain the origin of the LIGO-Virgo-KAGRA (LVK) observations lays between 1 and 100 M [44][45][46], which is precisely the region of parameter space of relevance for CMB anisotropy constraints on matter accretion onto PBHs depending on the underlying assumptions of the accretion process [27,28] (see, e.g., figure 2 of ref. [24] for a graphical representation).
Currently, the two most popular models of accretion mechanism onto a PBH assume spherical [27,47] or disk [28] accretion, with the latter leading to a significantly higher radiation luminosity than the former, hence to more stringent bounds. Similarly, also the type of ionization mechanism at play in the vicinity of the PBH has been shown to be able to affect the radiation luminosity by orders of magnitude [27]. However, in previous analysis the contribution from potentially present outflows has been neglected and it is therefore possible that their inclusion might alter the aforementioned conclusions. For instance, it is a priori unclear whether the dominant effect of the outflows is the reduction of the accretion rate via MF or the enhancement of the radiation luminosity due to the non-thermal emission of photons. The relative importance of these effects is very likely to depend not only on the PBH mass, but also on other assumptions regarding the geometry of the accretion, the ionization model and the characteristics of the outflows.
In this work we investigate these open questions by constructing an accretion model that accounts for both MF and non-thermal outflows, building on the analysis previously carried out in refs. [34,35]. We then study their impact on the thermal history of the universe and the CMB anisotropies in the context of both accretion geometries and ionization models, as well as for various PBH mass distributions. As a result, we find that the large theoretical uncertainties underlying the modelling of the outflows significantly enlarge those already present due to the other variables, leading to upper bounds on the PBH mass that can vary from approximately 0.01 to 100 M (assuming PBHs make all of the DM) depending on the chosen accretion model and outflow properties. This in turn bears important consequences for the LVK observability window. The work presented here shows that a more robust and complete modelling of PBH accretion is necessary to draw meaningful conclusions on PBHs being realistic DM candidates and progenitors of the LVK detected events. The paper is organized as follows. In section 2 we briefly introduce the generalities of PBH accretion and outline the current status of the field, while in section 3 we present the outflow model considered here in its details, adopting and extending the main results of [34,35]. In section 4 we succinctly explain how the accretion process affects the thermal history of the universe. In section 5 we discuss the numerical implementation of the accretion models, while in section 6 we calculate the resulting CMB anisotropy constraints for monochromatic and extended PBH mass functions, and discuss their implications for the LVK observed mass range. Finally, we conclude in section 7 with a summary of our analysis and closing remarks.

State of the art of PBH accretion
A PBH is moving supersonically in a homogeneously distributed gas accretes matter though an accretion column forming opposite to the direction of motion. In this scenario the accretion at the BH horizon is conveniently described by the Bondi-Littleton rate [48][49][50] where ρ ∞ is the mass density of the gas far away from the point mass, v rel is the relative velocity between the PBH and the gas, λ is the dimensionless accretion rate (which takes into account deviations from the idealised Bondi scenario due to the presence of e.g., pressure, viscosity, radiation feedback, MF, etc.), and r B is the Bondi radius, which characterizes the size of the BH sphere of influence and is defined as r B = GM PBH /v 2 rel . The consequent radiation luminosity of the system is typically parametrized as where is a dimensionless parameter that accounts for the radiation efficiency of the accretion process. 2 While the PBH mass is a free parameter of the model, quantities like v rel , λ and are determined by the properties of the environment close the PBH, such as, for instance, the temperature profile of the infalling gas, as well as by the geometry of the accretion or the details of the energy emission. Although there is no consensus on the exact form that these quantities should take, several ideas have been proposed so far in the literature (see e.g., refs. [27,28]). In order to set the stage for the accretion model described in section 3, where we focus more closely on the role of outflows, in this section we briefly review the state of the art of PBH accretion. Specifically, we focus on the role of the relative velocity between PBH and surrounding gas in section 2.1, on that of the gas ionization model in section 2.2 and on that of the geometry of the accretion in section 2.3, since these are the main aspects that affect the final accretion radiation luminosity in absence of outflows. The possible impact of other physical effects is discussed in section 2.4.

Accounting for PBH velocities
Constraints on the PBH abundance have been reported for different choices of the PBH velocities. Since the different assumptions can change the resulting constraints by orders of magnitude, we analyze and motivate here the approach employed in this work.
The first possibility is to consider that PBHs are at rest or moving at sub-sonic speed, hence the only relevant velocity is that of the surrounding medium, i.e., the speed of sound far away from the accretion region c 2 s,∞ = γv 2 B = γP ∞ /ρ ∞ , where γ = 5/3 is the polytropic equation of state index for an ideal monoatomic gas and v B = P ∞ /ρ ∞ is the Bondi velocity. Therefore, in the "v PBH = 0" scenario, where v PBH is the PBH proper velocity, we have v rel = c s,∞ . (2.3) However, under the assumption that the PBHs are the DM, PBHs are expected to be supersonic given that their velocity would effectively be the DM-baryon linear relative velocity v L , at least on large scales. 3 In this case the relative velocity reads Since velocities are stochastic variables, the radiation luminosity of a PBH population is in fact a velocity-averaged 4 luminosity L rad . So far two different averaging procedures have been presented in the literature [27]. In the first one, which we refer to as "approximated" average in the rest of the paper, the relative velocity is approximated as z being the redshift, and this result is used in equation (2.1) to obtain the radiation luminosity of the population. This approximation is well motivated in the cases where the radiation luminosity is proportional toṀ 2 PBH ∝ (c 2 s,∞ + v 2 PBH ) −3 , however in different scenarios its validity has to be assessed on a case-by-case basis. The approximation in equation (2.5) holds for all redshifts below z 10 4 and can therefore be safely applied for all times of interest for this work.
The second, more exact and general way of averaging over relative velocities is to consider that a transformation of the form c 2 s,∞ → c 2 s,∞ + v 2 PBH (which is the same as going from equation (2.3) to equation (2.4)) corresponds to a transformation of the gas temperature far away from the PBH of the form where m p is the proton mass andx e is the background free electron fraction. 5 In this case we express v rel in terms of T ∞ and average the resulting radiation luminosity over the distribution of the relative velocities to obtain the total radiation luminosity L rad . In the following we refer to this approach as "exact" average. We show in figure 1 the effect that these three choices for the relative velocity (v PBH = 0, approximated and exact average) have on the averaged total luminosity L rad . In the figure we consider the two types of ionization mechanisms described in section 2.2 as well as both the spherical and disk accretion scenarios detailed in section 2.3, for the representative case of M PBH = 10 2 M and f PBH = 1, where f PBH = ρ PBH /ρ DM is the fractional PBH energy density normalized to the DM abundance. 6 As can be clearly seen from the figure, the two methods to account for v rel when v PBH = 0 are overall almost equivalent, although differences (roughly of order two for the cases considered in the figure) can appear in particular for the disk accretion scenario, as in this case the radiation luminosity is not necessarily proportional toṀ 2 PBH . Moreover, we find that the exact approach is the most conservative. The figure further highlights the importance of including the PBH proper velocity, since the predictions of equation (2.3) always significantly deviate (even by orders of magnitude) from the exact approach. Therefore, for the sake of generality and being conservative, henceforth (unless stated otherwise) we always employ the exact averaging procedure to account for relative velocities, as done for instance in ref. [27].

Ionization models
In the scenarios of interest for this work, the radiation luminosity of an accreting PBH mostly comes from free-free radiation originated in the region close to the Schwarzschild radius [27]. The energy of the photons emitted in that region depends on the plasma temperature close to and around the compact object, which in turn depends on the dynamics of the infalling 5 This transformation follows from the fact that c 2 s,∞ = γT∞(1+xe)/mp, and hence T∞ = mpc 2 s,∞ /(1+xe)γ. 6 Although not explicitly present in the definition of the radiation luminosity, fPBH affects its value since L rad depends on the background temperature of the accreted matter, which is in turn sensitive to the total amount of injected energy and hence to fPBH, as discussed in section 4. We will therefore explicitly report the employed value of fPBH in all figures presented in the manuscript.
gas. Therefore, an accurate estimate of the temperature profile of the accreting region is necessary to determine the PBH radiation luminosity.
Before recombination, when the medium is already ionized, the temperature increases adiabatically the closer the gas gets to the BH. After recombination, when the medium far from the compact object becomes neutral, two scenario are possible, depending on the relative importance of two potentially competing effects: medium compression, which increases the temperature adiabatically, and ionization, which reduces the temperature proportionally to the energy lost in the ionization process and hence to the amount of initial neutral gas. In the first scenario the temperature of the infalling gas increases until it reaches T ion ∼ 10 4 K, then it stops increasing until all the gas gets ionized (since the balance between compression and ionization is maintained), and finally it resumes increasing. In this scenario the neutral gas is ionized through collisions with the free electrons and it is therefore referred to as "collisional ionisation". On the other hand, in the second scenario the radiation produced near the Schwarzschild radius is intense enough to directly photo-ionize the infalling gas before it reaches the ionization temperature T ion . Therefore in this scenario, which is conveniently called "photo-ionization", the temperature increases monotonically. In both scenarios radiative cooling is not expected to be relevant on scales close to the PBH.
Given that determining a realistic profile of the temperature around the PBH is very challenging and would require dedicated numerical simulations, we follow the simplifying approach of ref. [27]: we consider these two limiting scenarios of purely collisional or photoionization and report our results for both of them, knowing that they would bracket the "realistic" results.

Geometry of the accretion mechanism
The geometry of the accretion plays a crucial role in determining both the dimensionless accretion rate λ and the radiation efficiency . As in the previous section, also in this case we treat the spherical and disk accretion scenarios as the limits that encompass the "true" results.
In the spherical accretion scenario, most of the calculations can be performed semianalytically (see, e.g., ref. [27] for a recent treatment of the problem). The accretion efficiency can be expressed as a function of the Compton drag and cooling rate by CMB photons, finding limiting cases where these rates are either very efficient (isothermal accretion with λ 1.1) or negligible (adiabatic accretion with λ 0.1). It is also possible to study the evolution of the dimensionless accretion rate depending on the mass of the PBH and on the epoch of the universe, as can be seen in e.g., figure 4 of ref. [27]. Similarly, the radiation efficiency reduces to a function of the mass accretion rate and the temperature of the gas close to the PBH, which in turn depends on the ionization mechanism as discussed in the previous section. Typical values for the radiation efficiency in the spherical case are displayed in the left panel of figure 2 for difference choices of the PBH mass and both ionization models.
The other accretion mechanism we consider assumes the formation of a disk around the PBH [28]. Given that PBHs are expected to form binary systems already at early times [53], it seems natural to expect that the gas around the PBH experiences tidal forces that lead to the formation of a disk. Unfortunately, the details of disk accretion are very difficult to determine semi-analytically and only educated estimates have been used so far in the literature to determine quantities such as λ or . Following ref. [28], we fix as (arbitrary) benchmark value for the dimensionless accretion rate λ = 0.01, which has been argued to roughly take into account the role of viscosity and outflows close to the BH. 7 Similarly, we employ a parametric form of the radiation efficiency provided by ref. [54] (see in particular equation (11) and table 1 therein). However, there is an order-of-magnitude uncertainty in the estimate of the energy in electrons, characterized as the fraction of energy present in ions, which we assume here to have the representative value of δ = 0.1, as done in ref. [28]. The redshift evolution of in the disk accretion scenario for these benchmark choices is displayed in the right panel of figure 2. Note that these values are up to six orders of magnitude higher than those typical of the spherical accretion scenario. The figure also highlights that, while in the spherical accretion scenario ∝Ṁ PBH /L edd ∝ M PBH , in the disk accretion 6 is a parameter fitted to match numerical simulations of ref. [54], and hence the radiation efficiency scales differently as a function of the PBH mass in the two cases.
In summary, the unknown geometry of the accretion introduces an enormous source of uncertainty -potentially of orders of magnitude -in assessing the final radiation luminosity. This uncertainty is to be added on top of the uncertainty due to the ionization processes, although so far the disk accretion setup has been analyzed only in the context of collisional ionization. Due to these not completely understood physical processes, the constraints on the abundance of PBHs can vary substantially depending on the underlying choices made to describe the accretion process. Additional effects (discussed in detail in section 2.4) only further add to this uncertainty budget.

Additional effects
Besides these main factors, the accretion process can potentially be affected by several other contributions that have been considered in the literature. In this section we briefly outline a number of them, highlighting the regimes where they become relevant or can be neglected.
• The role of the Hubble expansion. This effect acts on the accretion process as a form of viscosity and introduces an additional redshift dependence of the background energy density and sound speed [55]. However, these contributions become important only for M PBH 10 4 M , threshold above which the quasi-steady flow approximation breaks down and so does the underlying mathematical setup. Therefore, in particular since in our analysis we focus on the LVK mass range, we do not consider masses above 10 4 M .
• The gravitational effect of the gas surrounding the PBH. The effect of self-gravity of the accreting gas can be neglected as long as the PBH mass is much larger than that of the surrounding medium, which is always the case for M PBH 10 5 M [55] (as for these masses the PBH can only attract a relatively low amount of matter) and therefore also for the range of PBH masses considered in this work. The same kind of consideration applies also for the gravitational impact of the background matter on the PBH (which could in principle affect the motion of the PBH, effect known as dynamical friction) as long as one considers PBH masses below 10 4 M [30,40].
• The role of super-Eddington accretion. Another possibility to consider for very massive PBHs is for them to accrete at a super-Eddington rate (see, e.g., refs. [56][57][58] and references therein for recent reviews). The Eddington accretion rate is the limit at which the gravitational force of the BH and the radiation pressure compensate (under the assumption of spherical symmetry), defining therefore the threshold for the maximum luminosity of the object. However, in practice, this limit can be exceeded for instance within proto-and massive galaxies at redshifts z 20, leading to periods of enhanced but variable luminosity. For spherical accretion, the transition to super-Eddington accretion is believed to happen for PBH masses roughly above 10 4 M , although in more realistic situations this value would strongly depend on the characteristics of the surrounding environment and on the geometry of the accretion (see, e.g., ref. [59]). Nevertheless, since this is a mainly low-redshift effect, our CMB bounds are largely unaffected by this source of uncertainty.
• The role of large-scale structure formation and DM halos. Once large-scale structures start to form at redshifts O(10), the non-trivial dynamics of PBHs and their interaction with the medium requires numerical simulations to be properly assessed [60][61][62][63], even if semi-analytical arguments can be made regarding how the accretion proceeds (see, e.g., ref. [64,65]). Furthermore, in the scenario where PBHs represent just a sub-dominant component of the DM and the rest is described by some new beyond-the-Standard-Model particle, DM halos are expected to form around the PBHs [66]. In this case, the formation of DM halos and of large-scale structures more in general can boost both the accretion of material [26,29,55] and the PBH proper velocities [28]. While the former effect enhances the radiation luminosity [30], the latter reduces it [65]. However, since these competing effects are strongly reliant on the assumptions made to describe the relevant late-time physics, we will neglect them in the current analysis for sake of clarity and leave their inclusion for future work.
• PBH mutual interactions. In the absence of large primordial non-Gaussianity, PBHs do not cluster over cosmic times [67] (see e.g., ref. [68] for further discussions) and PBH evolution can be effectively considered in isolation for masses below 10 4 M [27]. Furthermore, although early works suggested the existence of a large number of PBH binaries forming at early times [69], these are in contradiction with more recent numerical simulations [70,71]. Hence, PBH mutual interactions can be safely neglected in the mass range considered in this work.
• The role of the PBH spin. The radiation emitted by fastly rotating PBHs is enhanced by a factor of order ten percent with respect to the non-rotating case [72]. However, the majority of Kerr PBHs are realistically expected to have a dimensionless Kerr parameter much lower than unity [73,74], especially in the scenarios where PBHs form from the collapse of large primordial density fluctuations (see, e.g., refs. [75][76][77] for possible exceptions). Hence this effect can be neglected for the purpose of this work.
• The effect of inhomogeneuos energy injections. Recently, it has also been shown that the impact of PBH accretion on the photon bath can be very sensitive to the inhomogeneous spatial distribution of PBHs [78]. Nevertheless, since further important theoretical and numerical steps will be necessary to fully understand this effect, as also pointed out in the reference, we leave its inclusion in the analysis for future work.

The role of outflows
Another key physical effect that has been neglected in many of the previous analyses is the presence of outflows and their effect on large scales, where hereafter in this context "large" implies approximately the PBH scale of influence r acc 2GM PBH /v 2 PBH = 2r B . Recent analytical estimates and numerical simulations [34,35] suggest that even relatively weak jets and winds, i.e., outflows with different degree of collimation, might be able to sweep away some of the material around the PBH, thus reducing the flow of gas entering the BPH sphere of influence. Therefore, if present and carrying as little power as 10 −6Ṁ PBH [34], outflows could reduce the available material for accretion through this mechanical feedback effect. Despite the fact that the magnetic field in the primordial gas and the spin of the PBH are likely to be small, even pure hydrodynamical effects could drive outflows in a non-rotating black hole under small anisotropies of the accretion structure [34,[79][80][81]. Taking this and the uncertainties on accretion physics at very high z into account, the role of mechanical feedback deserves to be considered.
Furthermore, if outflows are produced during the accretion process, they can also potentially accelerate particles that can radiate and inject an additional amount of energy into the thermal bath. Such radiation from the outflows can compete with that produced in the accretion process itself, however, only for outflows that are relatively powerful. Therefore, outflows with a rather low power in terms ofṀ PBH could still produce mechanical feedback but generate negligible amounts of radiation. These two phenomena, analyzed separately in sections 3.1 and 3.2, are in fact competing effects: while decreasing the accretion rate intrinsically reduces the PBH radiation luminosity, such loss of luminosity is potentially compensated by the energy deposited into the medium by non-thermal particles accelerated by the outflows.
In this work we analyze the impact outflows potentially have on the two different accretion geometries presented in section 2.3. As explained before, these models have to be understood as limiting cases, presented to illustrate where the real physical scenario potentially lays in the PBH parameter space. Therefore, even if strictly speaking in the perfectly spherically symmetric scenario we cannot have outflows, it is still useful to consider that case as the benchmark reference for the scenario where some very low amount of angular momentum is present, which is expected in a cosmological scenario with small inhomogeneities, as outlined in ref. [34].

Impact on the accretion rate
As discussed in ref. [34], many different mechanisms promote the formation of outflows as, for instance, the existence of a magnetic field and/or the presence of some accreted gas with an excess of thermal energy. A (small) magnetic field can be generated by the Biermann battery mechanism, in which case the magnetorotational instability could enhance the magnetization and the turbulence, creating the required conditions for outflow formation [82].
In the presence of magnetic fields, outflows and jets can be fed by the rotational energy of the PBH or by the inner regions of an accretion disk [83,84], even when the accretion is quasi-spherical [85][86][87]. On the other hand, the production of an outflow from an excess of thermal energy in the gas is related to the formation of a thick and radiatively inefficient accreting disk, which forms winds through thermal pressure gradients. These outflows are expected to be broad and largely independent of the details of the magnetic field [32].
As demonstrated in ref. [34], once launched, outflows can overcome the ram pressure of the infalling gas and escape the PBH sphere of influence even when they are relatively weak and non-relativistic. Once outflows reach large scales, they compress and heat the medium by depositing energy and momentum into the gas. This phenomenon, called mechanical feedback (MF), is potentially responsible for a decrease in the accretion rate because the material is effectively heated and swept away from the PBH, implying that the amount of material that is available at scales of order r acc is reduced with respect to the standard case presented in section 2.
In order to account for this new large-scale effect, we modify the PBH accretion rate in equation (2.1) by introducing a new fractional rescaling f LS , so to obtaiṅ This new parameter specifically accounts for physical effects happening at large scales. Therefore, in the case with no MF, we have by definition that f LS ≡ 1 and we recover the cases presented in section 2. Conversely, when MF is present, f LS is expected to be lower than unity with a value that depends on the outflow ejected material velocity v o , the angle between the outflow and the PBH velocity θ and outflow half-opening angle χ (even if this latter dependence is weaker than the first two) [35]. Numerical simulations in 3D [35] suggest that a conservative choice for the fractional rescaling benchmark value (after averaging over all the possible values of the outflow orientation θ) is f LS = 0.1. Despite implying an order-of-magnitude reduction of the accretion rate, this choice is still cautious since phenomena such as the potential fast variation of the outflow orientation may further reduce the value of f LS and hence of the accretion rate (as discussed in [35]). Moreover, once PBHs become trans-or sub-sonic, the value of the fractional rescaling may be further reduced since the system becomes more axi-symmetric, and the accumulated internal energy and thus pressure within the outflow-medium interaction region can diminish accretion further (see, e.g., [34,88]). The in-depth analysis of such additional phenomena is left for future work.

Impact on the total luminosity
The presence of outflows might not only lead to MF, but also to the acceleration of nonthermal particles which later deposit energy and momentum into the medium. Therefore the total luminosity of the PBH receives contributions both from the radiation luminosity and the "non-thermal" luminosity induced by the outflow, i.e., where L rad is the same of equation (2.2) and the non-thermal contribution is conveniently parametrized by an additional efficiency parameter non-th as An upper-limit on non-th can be derived assuming that most of the energy released by the accretion gets released as outflows, which upon interaction with the environment (interaction by which MF is realized) converts most of its energy into non-thermal accelerated particles. In such a scenario, non-th would be ultimately determined by the relation between the dominant radiation channel, likely inverse Compton (IC) scattering off CMB photons, non-radiative losses via escape and adiabatic losses. Given the dense CMB photon field at z ∼ 10 2 − 10 3 , IC could efficiently cool relativistic electrons, with most of the radiation peaking somewhere within the range from soft X-rays to gamma rays. Such a scenario, which would imply non-th ∼ 0.1, is however highly optimistic, with more reasonable values for non-th 10 −2 (though values 10 −2 are also possible). A similar scenario has been explored in ref. [89] for the case of super-massive BH jets interacting with the medium at high z.
Due to the lack of quantitative predictions based on first principles, a robust (and not highly model-dependent) estimate of non-th has not been obtained yet. Therefore, in order to parametrize its impact on the cosmic medium, we consider three different benchmark values for non-th , namely non-th = 10 −2 , 10 −4 , 10 −6 . As a reference, the intermediate case, non-th = 10 −4 , may be interpreted as a scenario where 10% of the accretion energy goes to the outflows, of which 10% goes to non-thermal particles, which in turn can radiate through IC about 10% of their energy, 10% of which in the form of photons with suitable energies (see ref. [86] for a discussion in the context of isolated black holes in our galaxy). These benchmark values are compared to the radiation efficiencies obtained for the scenarios in absence of outflows in figure 2. As it can be observed in the figure, the chosen non-thermal emission efficiencies are of the same order of magnitude (depending on the PBH mass) of the radiation efficiencies obtained in the represented cases, hinting to a tight interplay between the two terms contributing to equation (3.2).
Finally, we note that in this case the gas can be easily ionized through compression by the shock generated by the outflow on large scales. Hence, since the gas cooling time is not small enough for the gas to recombine, the photo-ionization model provides a more realistic description of the physical picture. Nevertheless, we will also consider the collisional ionization model for sake of completeness.

Combined effect
The impact of MF and non-thermal emissions on the total luminosity is shown in figure 3 for both accretion geometries (spherical on the left, disk on the right) and ionization models. For sake of succinctness (but with no loss of generality), here we assume f PBH = 10 −3 and non−th = 10 −4 . Furthermore, we also show the results only for the representative case of z = 600, which has been shown to be the time around which the CMB anisotropies are most responsive to (almost constant in time) energy injections such as the ones considered here [90]. A more general redshift dependence can be qualitatively inferred from figure 2.
In figure 3, the various predictions in absence of outflows are shown in blue and match the expected linear dependence of the radiation luminosity on M PBH (with different powers depending on the geometry as explained in section 2.3) over almost the full PBH mass range. These curves are then suppressed by orders of magnitude when including the role of MF 8 , as represented by the magenta lines. However, this reduction of luminosity can be partially compensated, in particular in the low mass range, when non-thermal emissions are taken into account (red lines).
This implies that when including the effect of outflows there are three possible outcomes.
1. First of all, if non-thermal emissions are very efficient, their contribution to the total luminosity always dominates over the radiation luminosity and results in an enhanced energy injection with respect to the case in absence of outflows (unless very low f LS values, say 1, are considered). This scenario is represented, for instance, in the left panel of figure 3 in the collisional ionization case.
2. On the other hand, the opposite is true if non-thermal emissions are relatively inefficient, as observable in the disk accretion case represented in the right panel of the figure.
3. Finally, as can be seen from the left panel of the figure for the photoionization case, there might be a range of PBH masses where the total luminosity is enhanced with respect to the case in absence of outflows due to non-thermal emissions (low mass range) and regions where this is instead suppressed because of MF (high mass range). 8 That is, since L rad ∝Ṁ 2 PBH ∝ f 2 LS in the spherical accretion case and L rad ∝Ṁ 1+a PBH ∝ f 1+a LS in the disk accretion case. By assuming fLS = 0.1 (as argued in the previous section) one obtains a suppression of a factor 100 and 40, respectively.
As discussed in section 6, the balance between these regimes dictates how the final CMB constraints are affected by the presence of outflows.

Impact on thermal history and CMB anisotropies
The last ingredient to factor in our analysis is the description of how the emitted radiation impacts the thermal history of the universe. A population of PBHs injects energy into the cosmic medium at a rate whereρ PBH is the PBH background energy density, f PBH =ρ PBH /ρ cdm is its fractional form with respect to the DM abundance, and L tot is the velocity-averaged total luminosity. However, the amount of energy that is actually deposited into the thermal bath can in principle be different, for instance depending on the transparency and energy density of the medium or on the type of emitted particles [91,92]. Furthermore, different deposition channels have to be accounted for, since the additional energy might not only heat up the medium but also ionize and excite it if the injection happens after recombination [93][94][95].
Following [96][97][98], we express the deposited energy rate in a channel c as where f eff and χ c are the deposition efficiency and deposition fraction for the specific channel, respectively. 9 In this work we parameterize f eff as in section IV of ref. [27], while for the calculation of the different deposition fractions we rely on table V of ref. [95]. Once the relation between injected and deposited energy is set, one can analyze how the different accretion scenarios affect the reionization fraction x e , which is ultimately the quantity that determines how much the CMB anisotropy power spectra are affected by the non-standard energy injection. In fact, the main effects of depositing energy into the cosmic medium are twofold: (i) delaying the time of last scattering, which induces a slight shift in the position of the peaks in the CMB temperature and polarization power spectra, and (ii) effectively increasing the ionization optical depth since the reionization of the universe starts at an earlier time. In figures 4 and 5 we highlight this behavior by comparing the effect on the free electron fraction x e and the CMB anisotropy power spectra (computed as described in section 5) of the different accretion models considered in this work, with and without outflows, for three representative choices of PBH mass.
From the figures, several conclusions can be drawn. On the left panels of figure 4, which represent the spherical accretion case, we compare the cases with and without outflows and we find that for low PBH masses (upper panel) the contribution from the outflows always dominates over the radiation luminosity. In fact, once MF is accounted for, the mass accretion rate is suppressed by a factor 10 (assuming the fiducial values introduced in section 3.1), and the corresponding radiation efficiency becomes significantly lower than any of the cases with a non-zero non−th . Interestingly, this is true for both ionization models discussed in section 2.2, so that no difference is to be seen between the collisional (dashed lines) and photoionization (solid lines) cases when including outflows in the low mass range. The situation  . Role of the outflows in the impact on the free electron fraction of the accretion process (shown are the free electron fraction -lower panels -and its relative difference with respect to the ΛCDM prediction -upper panels). From top to bottom we compare three representative choices of the PBH mass (same as in figure 2), while from left to right we compare the spherical and disk accretion scenarios (for the latter we assumed λ = 0.01 and δ = 0.1). In the bottom panels the ΛCDM prediction is shown in black. The irregular shape of the collisional ionization curve in the top left panel is due to numerical noise.   slightly changes, however, for PBH masses of the order of 10 2 M (middle panel) and nonthermal efficiencies lower than 10 −4 , where the radiation luminosity in the photo-ionization case becomes comparable to that of the non-thermal outflows after the inclusion of MF and differences between the various ionization processes start to emerge. This behavior is further enhanced for larger PBH masses (bottom panel), although the non-thermal luminosity still dominates for non−th ∼ 10 −2 .
On the other hand, because of the higher radiation efficiencies reached in the disk accretion scenario (see figure 2), already at relatively low PBH masses (upper panel) the intrinsic radiation efficiency of the accretion is comparable to the non-thermal contribution once MF is taken into account for non−th ∼ 10 −6 . For higher PBH masses (middle and bottom panels) the role of the non-thermal emission becomes negligible and the curves for non−th < 10 −4 perfectly overlap, with only a minor difference in the non−th = 10 −2 case.
In summary, we expect the contribution from the outflows to be largely determining the luminosity of the system for low masses in the spherical accretion scenario, greatly reducing (if not erasing) the differences between the various ionization models in this regime, while for larger PBH masses the interplay between radiation luminosity, MF and non-thermal emissions starts to become apparent. This differs from the disk accretion case where the intrinsic radiation luminosity of the system is always relevant, and eventually dominates already for PBH masses of the order of 10 2 M .
In figure 5 we also show the impact that these modified thermal histories have on the CMB temperature and polarization anisotropy power spectra. For sake of succinctness, here we focus only on the M PBH = 10 2 M case, which we believe already suffices to convey the main message of the section. In particular, for both accretion geometries, we notice the same level of interplay between radiation efficiency, MF and non-thermal effects as in the corresponding panel of figure 4. Interestingly, the contribution of the outflows (visible for instance in the non−th = 10 −2 case of the left panels) resembles very closely the smooth behavior already observed in ref. [28] for the disk accretion scenario (see figure 3 of the reference), while the purely spherical accretion case imprints a more oscillatory behavior in the residuals (see, e.g., figure 13 of ref. [27]). This might mean, for instance, that a potential observation of a similar signal might be unable to disentangle a spherical accretion scenario dominated by outflows from the disk accretion scenario, regardless of the impact of outflows.

Numerical setup
All the relevant cosmological quantities for this work, such as the CMB power spectra, are computed employing the latest version of the Boltzmann solver CLASS [101,102]. In particular, we make use of the energy injection treatment discussed in ref. [96], which is largely based on the ExoCLASS extension of CLASS [99]. Recombination is solved using the HYREC [103,104] implementation of CLASS, while the energy injection and deposition process is computed using the prescriptions detailed in section 4.
In terms of PBH accretion, we further improve upon the implementation of ref. [96] (and [99]) by including also the photo-ionization option for the spherical accretion case as well as the calculation for the average of the PBH luminosity over the PBH proper velocities, which we employ as default option for our analysis (see section 2.1 for further details). The impact of the outflows is implemented as an extension of the spherical and disk accretion cases, where we simply modify the definitions of the accretion rateṀ PBH and of the total luminosity L tot , as illustrated by equations (3.1) and (3.2), respectively. When considering disk accretion, we fix λ = 0.01 and δ = 0.1, the benchmark values suggested in ref. [28].
The constraints on the cosmological parameters discussed in section 6 for the different accretion models are obtained using the parameter inference code MontePython [105,106]. As commonly done in the literature, we perform a number of Markov Chain Monte Carlo (MCMC) scans of the parameter space for fixed values of the PBH mass, which reduces the problem to a 6+1 extension of the ΛCDM model with where ω b and ω cdm are the baryon and DM physical energy densities, respectively, h is the dimensionless Hubble rate, A s and n s are the amplitude and scalar spectral index of the primordial power spectrum, and τ reio is the reionization optical depth. 10 In the scenarios of interest, PBHs are described by the abundance parameter f PBH =ρ PBH /ρ cdm , i.e., by the fraction of DM in form of PBHs. We consider temperature, polarization and lensing information from the Planck 2018 mission [1] (explicitly, we use the high-TTTEEE, low-EE, low-TT and lensing likelihoods) and determine the MCMCs to be converged with the Gelman-Rubin criterium |R − 1| < 0.02 [107].

CMB constraints on PBH abundance
In this section we derive the corresponding CMB anisotropy constraints on the PBH abundance for both accretion geometries and ionization models described in section 2 and focusing on the impact of outflows for different choices of the non-thermal emission efficiency as explained in section 3. To ease the comparison with the literature, first of all we derive the aforementioned constraints in the context of a monochromatic PBH mass distribution in section 6.1. Then, for sake of generality, we also recast our constraints in terms of a popular choice of extended mass distribution in section 6.2. On the basis of these results, in section 6.3 we also comment on their implications for the LVK window.

Constraints for monochromatic mass distributions
We begin by assuming that PBHs have a monochromatic mass distribution (MMD), i.e., all PBHs have the same mass. In full generality, a PBH mass distribution is described by the fractional abundance function [108] df where dΦ PBH /dM PBH describes the shape of the PBH mass distribution, it is normalized to unity by construction and in the monochromatic case it reads as In this work we choose M PBH such that it scans the mass range [10 −2 , 10 4 ] M , computing for each PBH mass the respective upper limit on the PBH abundance. We show in figure 6 the corresponding 95% CL upper bounds on the fractional PBH abundance f PBH for different accretion mechanisms, ionization choices and outflows scenarios. In particular, in each of the two subplots, which display the bounds for the spherical accretion case on the left and the disk accretion case on the right, we report the upper bounds on f PBH for the two aforementioned ionization models (photo-ionization reported as solid lines and collisional ionization as dashed lines) as well as for different choices of the parameters describing the outflow modelling. As in the previous sections, we have the scenario without any outflow contribution in blue, while the contours including MF with the benchmarking value of f LS = 0.1 and non-thermal emissions are reported in magenta, red and orange for non−th = 10 −6,−4,−2 , respectively.
Focusing first on the left panel of figure 6, i.e., on spherical accretion, we observe that the constraints obtained in absence of outflows are approximately one order of magnitude more stringent than those reported in ref. [27] for both ionization models. We attribute this discrepancy to the difference in methodology 11 and Planck data release employed here compared to the one in ref. [27]. The bounds for the collisional ionization case are in fact broadly  Figure 6. Impact of outflows on the CMB constraints on the fractional PBH abundance for the two main types of accretion geometry (spherical on the left and disk on the right) and ionization models (photo-ionization in solid and collisional ionization in dashed) assuming a MMD. The blue contours represent the scenarios without outflows and update the bounds first derived in [27] (left) and [28] (right), while the magenta, red and organge lines assume MF with f LS = 0.1 and non-thermal emissions with non−th = 10 −6,−4,−2 , respectively. consistent with those presented in ref. [30], where a more similar analysis was conducted. In terms of impact of the outflows, for low non-thermal efficiencies the role of MF alone is dominant, leading to an overall suppression of the constraints by approximately 1 − 2 orders of magnitude, as explained in section 3.3. For intermediate values of non−th the interplay between MF and non-thermal emissions becomes more balanced, with the latter most predominately contributing to the total luminosity until PBH masses of the order of 10 3 M , above which the radiation luminosity dominates despite the presence of MF. Nevertheless, also in this scenario we observe a suppression of the bounds with respect to the case without outflows, with only a minor enhancement for masses below 10 M . The situation changes for very efficient non-thermal outflows, which dominate the accretion emission over the whole PBH mass range, thereby arasing the difference between photo-and collisional ionization (which perfectly overlap in the figure), and significantly strengthen the constraints down to PBH masses of the order of 0.1 M for f PBH = 1.
On the other hand, since in the context of disk accretion the intrinsic radiation luminosity is much higher than in the spherical accretion case, in the right panel of figure 6 we notice that even for non-thermal efficiencies as high as 10 −4 the only relevant impact of outflows is MF, which suppresses the constraints by about 1 − 2 orders of magnitude. For this reason, the contours for the non−th = 10 −4 and non−th = 10 −6 cases perfectly overlap in the right panel of the figure. Only for non−th of the order of 10 −2 non-thermal emissions start to play a significant role, enhancing the constraints for PBH masses below 10 M but still leading to a suppression thereof for larger PBH masses.
Focusing on the shapes of the constraints, one expects the behaviour of the bounds to follow the dependence of the deposited energy on the free parameters of the model, i.e., f PBH and M PBH , which roughly boils down to dE/dtdV | dep ∝ f PBH M PBH λ( non−th + M 2 PBH λ). In the case of spherical accretion the values of λ at z ∼ O(10 − 10 3 ) increase as a function of the mass (see figure 4 of ref. [27]), which explains why the relation between PBH abundance and mass is not linear, with the bounds becoming comparatively less stringent the lower the PBH mass (a similar behaviour is also to be observed in figure 3). This relation, and in particular the different dependence on the PBH mass, also explains the difference in slope of the f PBH − M PBH upper bounds between the radiation and non-thermal emission dominated regimes. In the disk accretion scenario the value of λ is fixed to a fiducial value, and therefore the upper bounds approach a linear dependence more closely than in the spherical case. Nevertheless, as it becomes clear in particular in the non−th = 10 −4,−6 cases, the parametric form of used in ref. [54] still introduces a non-trivial dependence of the constraints on the PBH mass.

Constraints for extended mass distributions
Even if for practical purposes it is more convenient to obtain constraints for MMDs, in reality it is well known that PBH populations would likely have an extended mass distribution (EMD). Different PBH formation mechanisms are generally responsible for different EMDs, although we can identify two popular benchmark classes of EMD: power-law and lognormal. Broadly speaking, the first class is typically associated with the collapse of large density perturbations or cosmic strings, see e.g., refs. [109,110], while the latter is connected to the presence of large peaks in the primordial power spectrum, see e.g., refs. [111][112][113].
Recent works, as for instance refs. [114][115][116][117][118][119][120][121][122][123][124][125][126], have shown significant interest in models where PBHs are generated from peaks in the primordial curvature power spectrum. Such effort has been mainly driven by the intrinsic connection between the primordial power spectrum and inflationary dynamics, which is known to be compatible with the single-field slow-roll scenario only during 6/8 e-folds of the at least 50 required to solve the horizon and flatness problems [127][128][129]. Therefore establishing robust bounds on the PBH abundance, or even just their presence, is a key ingredient in reconstructing the primordial power spectrum at scales of order O(10 5 − 10 15 ) Mpc −1 currently not accessible by other cosmological observables [130][131][132][133][134][135][136][137][138][139][140]. For this reason, here we focus on deriving abundance constraints for lognormal EMDs characterized by two parameters, the mean µ and standard deviation σ. More realistic mass distribution models would require a greater number of parameters to describe the EMD [141], and therefore in this work we choose the simplest extension of the MMD case. The same procedure can be implemented also in the case of power-law EMDs, as done for instance in ref. [23].
Interpreting abundance constraints obtained for MMDs as abundance constraints for EMDs is not straightforward [142], although several methods have been proposed to infer the latter from the former (see e.g., refs. [108,142,143]). In this work we follow the approach of ref. [108], based on the concept of equivalent mass. In the case of interest, if PBHs have an EMD, equation (4.1) reads as showing that it is possible for different EMDs to inject the same amount of energy into the cosmic medium, making it possible for PBH populations in those models to be equally abundant. At a practical level, it is always possible to associate the impact on a cosmological observable of a given EMD, and therefore its abundance f EMD PBH , to the same-magnitude effect generated by a MMD with equivalent mass M eq , i.e., to its abundance f MMD PBH . Despite its  Figure 7. Same as in figure 6, but for a lognormal EMD. The figure shows the collisional and photoionization regimes for the spherical accretion case (left and center panels, respectively) and the disk accretion model (right panels). We report here all the cases in absence of MF and non-thermal emission (top panels), and with those effects included (middle and bottom panels) for different nonthermal efficiencies. The grey shaded area represents the region of parameter space that cannot be described appropriately by the theoretical model employed in this work. simplicity, this method provides an easy analytical tool to compute constraints for EMDs starting from the MMD analysis of section 6.1. The validity of this method has been already explicitly proved in ref. [23] for the CMB case or, for instance, in ref. [38] for other abundance constraints also based on accretion physics. In figure 7 we show the PBH abundance constraints for the lognormal EMD for all geometries, ionization models and non-thermal emission efficiencies already employed in the previous section. We report the new conversion formula used for all the cases at hand in appendix A. The gray shaded regions in the plots indicate the values of µ and σ for which the EMDs extend beyond 10 4 M and the theoretical models employed in this work break down as discussed in section 2.4 (see ref. [108] for further details). By comparing the cases with only radiative efficiency to the ones with also non-thermal emission, we note that as soon as the MF effect starts to play a role, the allowed parameter space increases significantly with respect to the case of radiative emission only. On the other hand, as soon as the non-thermal emission grows, the parameter space becomes very tightly constrained, as can be seen in the bottom panels. In other words, because of these two competing effects, the same EMD can be either ruled in or ruled out depending on the details of the accretion and emission models.

Theoretical uncertainties and implications for the LVK mass range
The analysis presented in this work clearly highlights what is one of the underlying issues with existing constraints on PBHs: the unknown magnitude of the theoretical error associated to the modelling of the accretion process. Reducing this uncertainty is crucial given that existing and future GW observatories are and will be sensitive in the O(10 − 10 4 ) Hz frequency range, i.e., the typical frequency range where the signal of coalescing compact objects with masses O(10 − 10 3 ) M is expected to be. An accurate estimation of the PBH abundance upper limits would in fact allow us to speculate more realistically about what amount of detected events can have a primordial origin.
Nevertheless, such an accurate description of the accretion process is still not available. Therefore, we present our findings in figure 8 as "uncertainty bands" around the true, unknown, upper bounds. These bands encapsulate the effects of the geometry and ionization models with (orange) and without (blue) the inclusion of outflows for a MMD. As shown in the figure, even before accounting for the outflow phenomenology, the width of the uncertainty band spans over two orders of magnitude, approximately between 1 and 100 M , for f PBH 1. This means that the intrinsic theoretical error in the modelling of the accretion makes these CMB anisotropy constraints unable to discern whether PBHs might be the progenitor of any observed BH merger, as they could either completely exclude or allow the LVK mass range. Once outflows are included in the modelling, the uncertainty band significantly enlarges both at large PBH masses, where MF dominates and the effect of non-thermal emission is negligible, and at low PBH masses, where the effect of non-thermal emission can be sizeable. Furthermore, the uncertainty due to the theoretical modelling in the case of EMDs can be easily seen when comparing the different panels of figure 7.
In this respect, we also note that the existing bounds on the PBH abundance deriving from the LVK estimate of the binary BH local merger rate are subject of an ongoing debate. In fact, semi-analytical estimates of the local merger rate R 0 of PBH binaries formed at early times suggest that it could be as high as R 0 ∼ 10 5 Gpc −3 yr −1 [11,69] for f PBH = 1, to be compared with the observational value of R 0 ∼ 20 Gpc −3 yr −1 [144]. On the other hand, recent (and more accurate) numerical simulations reconcile the PBH early binary local merger rate value with the observed one for f PBH = 1 [70,71].
The question of whether PBHs can be a significant component of the DM content of the universe and of the GW events detected by LVK is therefore still open and only a synergistic effort between different communities will help us find an answer to it. In fact, although in this paper we focused only on CMB constraints, they are not the only relevant ones in the LVK mass range. For instance, there are many complementary constraints coming from supernova lensing [145], dwarf and ultra-faint dwarf galaxy dynamics [146][147][148], Lymanα forest [149,150] or wide binaries survival [151], which might help us shed light on the complexities of PBH phenomenology and ultimately determine their abundance.

Conclusions
Despite its remarkable success in explaining numerous cosmological observables, the ΛCDM model cannot provide any insight on what the true nature of DM is, allowing a vast plethora of models to fit observations. For the sake of convenience, DM candidates are typically categorized in terms of their mass: already in this context, DM models span many orders of magnitude, ranging from ultralight axion-like particles to macroscopic compact objects like BHs. Hence finding novel ways to constrain DM properties becomes fundamental in order to establish its nature.
In this work we focus on one of these popular candidates, PBHs. In particular, we focus on PBHs with masses larger than 1 M , since they have the potential to be detected by existing and future GW observatories. This class of PBHs had been previously thought to be ruled out by existing LVK constraints on the local merger rate, although recent numerical simulations have shown that the number of PBH merging binaries had been overestimated by several orders of magnitude, reopening the possibility for PBHs to be a substantial component of the DM in the LVK mass range.
In this spirit, we turn our attention to other existing constraints on the PBH abundance in that same mass range, that is to those coming from the accretion of matter into a PBH. In fact, the emission of radiation following the accretion process affects the thermal history of the universe, by delaying recombination and anticipating reionization, and can therefore be constrained by CMB observations. However, large theoretical uncertainties underlay the modelling of accretion and this translates in large error bars on the final constraints. Examples of sources of such uncertainties are the geometry of the accretion as well as the ionization model determining the temperature profile close to the BH.
On top of these uncertainties, one aspect of the accretion physics that has not been considered extensively in the literature so far is the effect of outflows (winds and/or jets, depending on their degree of collimation) on the accretion. However, it has been shown both analytically and numerically that even relatively weak outflows can escape the BH sphere of influence and sweep away at least part of the cosmic medium around the PBH, thereby decreasing its accretion rate. At the same time, the very same outflows could also accelerate non-thermal particles and effectively enhance the luminosity of the BH. Therefore, these competing effects can significantly affect the total luminosity of the system and introduce an additional layer of uncertainty that needs to be taken into account when quoting cosmological constraints derived from PBH accretion.
In this work we attempt to model the largely unknown nature of these effects and to analyse how their balance affects the CMB constraints on PBH accretion. We do so for different choices of accretion geometry (spherical or disk) and ionization models (photo-and collisional ionization), as well as for both monochromatic and extend mass distributions. Our quantitative findings are shown in figures 6-7, which clarify that the final outcome heavily relies on the choice made in particular for the efficiency of the non-thermal emissions. This conclusion is graphically summarized in figure 8, where the cumulative uncertainty bands on the "true" bounds are shown with and without the inclusion of outflows.
Looking towards the future, several developments might have a significant impact on our results. On the one hand, should strongly compelling arguments be put forward in favour of one particular accretion geometry, it would significantly reduce the uncertainty on the PBH abundance upper bounds reported in figure 8. The same would be true also in the context of the ionization model. On the other hand, however, even more accurate simulations than the ones performed in [34,35] might find indications for more complex outflow dynamics requiring, for instance, also information on outflow orientation or the transition to transor sub-sonic regimes. Taking into account for the potential role of these unknowns would inevitably further widen the size of the aforementioned uncertainty band. Similarly, the disk accretion scenario considered here is only restricted to the collisional ionization case, and extending it to the photo-ionization model would extend the uncertainty region towards low PBH masses.
Overall, we conclude that the path towards a realistic estimate of whether PBHs can make up for a sizeable fraction of the DM and of the events observed at the LVK facilities at the same time (i.e., in the same mass range) is still very long, as both the cosmological constraints and the merger rate estimates are subject to very large theoretical uncertainties. While this should discourage premature claims of exclusion or detection, it should also be seen as a source of motivation for the improvements to come. A Converting abundance constraints from MMD to EMD Abundance constraint conversion formulas accounting for PBH effects on CMB anisotropies exist only for the spherical accretion case [36,108]. However, different accretion geometries or energy emission efficiencies require the development (and testing) of new conversion relations, which we present in this appendix. Furthermore, in this appendix we focus on lognormal EMDs, given the large interest on this specific EMD coming from the theoretical modelling of PBH formation. We note, however, that the procedure presented in this appendix is easily applicable also to other EMDs.