Multi-wavelength astronomical searches for primordial black holes

TL;DR

This paper analyzes radio and X-ray signals from primordial black holes in our galaxy to set new, stronger limits on their contribution to dark matter, improving previous bounds with advanced accretion physics modeling.

Contribution

It introduces an improved gas accretion model based on numerical simulations and relaxes the monochromatic mass assumption, refining constraints on primordial black hole dark matter.

Findings

Maximum relic density of primordial black holes is about 0.1% of dark matter.

New upper bounds are two orders of magnitude stronger than previous estimates.

Future multi-wavelength searches will better probe galactic black hole populations.

Abstract

If primordial black holes of $\mathcal{O}(1-100) M_{\odot}$ constitute a significant portion of the dark matter in the Universe, they should be very abundant in our Galaxy. We present here a detailed analysis of the radio and X-ray emission that these objects are expected to produce due to the accretion of gas from the interstellar medium. With respect to previous studies, we relax the assumption of a monochromatic mass function, and introduce an improved treatment of the physics of gas accretion onto isolated, moving compact objects, based on a set of state-of-the-art numerical simulations. By comparing our predictions with known radio and X-ray sources in the Galactic center region, we show that the maximum relic density of primordial black holes in the mass range of interest is $\sim 10^{-3}$ smaller than that of dark matter. The new upper bound is two orders of magnitude stronger with respect to previous results, based on a conservative phenomenological treatment of the accretion physics. We also provide a comprehensive critical discussion on the reliability of this bound, and on possible future developments in the field. We argue in particular that future multi-wavelength searches will soon start to probe the galactic population of astrophysical black holes.

Figures

Equations

$$v_{\mathrm{mean}}=v_{\mathrm{circ}}(R)=\sqrt{\left(GM(<R)/R\right)}.$$

$$\dot{M}=\lambda {\dot{M}}_{\mathrm{BHL}},$$

$${\dot{M}}_{\mathrm{BHL}}\equiv 4\pi (G{M}_{\mathrm{BH}}{)}^2\rho (v^2+c_{\mathrm{s}}^2{)}^{-3/2}$$

$$L_B=\eta \dot{M}{c}^2$$

$${\dot{M}}_{\mathrm{crit}}\equiv 0.01\cdot \frac{L_{\mathrm{Edd}}}{{\eta }_0{c}^2},$$

$$f(M|\mu ,\sigma )=\frac{1}{M\sqrt{(}2\pi )\sigma }\exp \left(-\frac{log(M/\mu {)}^2}{2{\sigma }^2}\right),$$

$$f(M|\gamma ,M_{\min },M_{\max })=\frac{\mathcal{N}}{M^{1-\gamma }}\Theta (M-M_{\min })\Theta (M_{\max }-M),$$

$$\mathcal{N}=\{\begin{array}{ll}\frac{1}{log(M_{\max }/M_{\min })},& \text{if}\gamma =0,\\ \frac{\gamma }{M_{\max }^{\gamma }-M_{\min }^{\gamma }},& \text{if}\gamma \ne 0.\end{array}$$

$$\gamma =-\frac{2w}{1+w}.$$

$$f(M)\propto \{\begin{array}{ll}{M}^{-1},& \text{if}{M}_{\min }<M<M_{\ast }\\ M^{-3/2},& \text{if}M>M_{\ast }\end{array}$$

$$f_{\mathrm{PBH}}\int g(M,\{p_j\})\frac{d\Phi }{dM}dM,$$

$$f_{\mathrm{PBH}}^{\mathrm{MMD}}g(M_{eq},\{p_j\})=f_{\mathrm{PBH}}^{\mathrm{EMD}}\int g(M,\{p_j\})\frac{d\Phi }{dM}dM,$$

$$\langle\lambda_{\mathrm{B}}\rangle=\frac{\rho_{\mathrm{in}}(GM_{\mathrm{BH}})^{2}}{(v^{2}_{\mathrm{in}}+c^{2}_{\mathrm{s,in}})^{3/2}}\frac{c^{3}_{\mathrm{s}}}{\rho(GM_{\mathrm{BH}})^{2}}=\frac{\rho_{\mathrm{in}}}{\rho}\Big{(}\frac{c_{\mathrm{s}}}{c_{\mathrm{s,in}}}\Big{)}^{3}\frac{1}{(1+v^{2}_{\mathrm{in}}/c^{2}_{\mathrm{s,in}})^{3/2}},$$

$$\dot{M}=⟨{\lambda }_{\mathrm{B}}⟩{\dot{M}}_B={\Delta }_{\rho }(2{\Delta }_T{)}^{-3/2}(1+{\mathcal{M}}_{\mathrm{in}}{)}^{-3/2}\frac{\pi \rho (G{M}_{\mathrm{BH}}{)}^2}{c_{\mathrm{s}}^3}{e}^{3/2},$$

$${\mathcal{M}}_{\mathrm{in}}=\frac{v_{\mathrm{in}}}{v_{\mathrm{BH}}}\frac{v_{\mathrm{BH}}}{c_{\mathrm{s}}}\frac{c_{\mathrm{s}}}{c_{\mathrm{s},\mathrm{in}}}=\frac{1}{{\Delta }_{\rho }}\mathcal{M}\frac{1}{(2{\Delta }_T{)}^{1/2}}.$$

$${\Delta }_{\rho }=\{\begin{array}{ll}\frac{1+{\mathcal{M}}^2+\sqrt{(1+{\mathcal{M}}^2{)}^2-8{\mathcal{M}}^2{\Delta }_T}}{4{\Delta }_T},& \mathrm{if}\mathcal{M}<{\mathcal{M}}_D,\\ \frac{1+{\mathcal{M}}^2}{4{\Delta }_T},& \mathrm{if}{\mathcal{M}}_D\le \mathcal{M}\le {\mathcal{M}}_R,\\ \frac{1+{\mathcal{M}}^2-\sqrt{(1+{\mathcal{M}}^2{)}^2}-8{\mathcal{M}}^2{\Delta }_T}{4{\Delta }_T},& \mathrm{if}\mathcal{M}>{\mathcal{M}}_R.\end{array}$$

$$N_S={\int }_{{\varphi }_{\mathrm{thres}}}^{\infty }\frac{dN}{d\varphi }d\varphi ={\int }_0^{\infty }g(M,\{p_j\})\frac{dN}{dM}dM.$$

$$g(M,\{p_j\})={\int }_{{\varphi }_{\mathrm{thres}}}^{\infty }p(h(\varphi ,M))\left|\frac{\partial h(\varphi ,M)}{\partial \varphi }\right|d\varphi .$$

$$g(M,\{p_j\})=\{\begin{array}{ll}0,& \mathrm{if}M\le M_{\min },\\ \mathrm{erf}\left(\frac{h_{\mathrm{thres}}}{\sqrt{2}a}\right)-\sqrt{\frac{2}{\pi }}\frac{h_{\mathrm{thres}}}{a}\exp \left(-\frac{-h_{\mathrm{thres}}^2}{2a^2}\right),& \mathrm{if}M>M_{\min },\end{array}$$

$$h({\varphi }_R,M)={\left({\varphi }_R^{-1.45/3}X{M}^{3.88/3}-c_s^2\right)}^{1/2}\mathrm{and}h({\varphi }_X,M)={\left({\varphi }_X^{-1/3}YM-c_s^2\right)}^{1/2}$$

Full text

1 Introduction

The idea that primordial black holes (PBHs) exist in the Universe [1, 2, 3] has been widely discussed in the past few decades [4]. These hypothetical objects could have formed before big-bang nucleosynthesis out of the collapse of small-scale large-amplitude density fluctuations originated during inflation, or through a variety of other mechanisms including phase transitions, topological defects such as cosmic strings and domain walls, condensate fragmentation, and bubble nucleation (see e.g. Refs. [5, 6]).

PBHs with mass between $\sim 10^{-18}$ and $\sim 10^{6}$ solar masses have been invoked as possible dark matter candidates [2, 7]. One of the most interesting mass window appears to be the one associated with the massive black-hole-binary merger events reported by the LIGO and VIRGO collaborations, centered on $\mathcal{O}(10)$ solar masses. These events have indeed significantly revived the interest in this topic [8, 9, 10], encouraged a re-evaluation of the existing bounds [11, 4, 12], and triggered an extended discussion on the prospects for astronomical and cosmological signatures [13, 14].

In a previous paper [15] we have shown that the hypothesis that PBHs comprise all the DM in the universe can be constrained by exploiting Galactic astronomical data, in the radio and X-ray bands: The idea is that a small fraction of these objects are expected to accrete interstellar gas in a significant way, especially in the inner Galactic bulge, and therefore show up as radio and X-ray sources in the sky. By comparison with current source catalogs we could place a conservative bound on their abundance in the $10$ - $100$ $M_{\odot}$ mass window.

In view of the upcoming observing runs of advanced gravitational waves interferometers, and in preparation for the exquisite radio observations that will become available thanks to the SKA project – an array of detectors that will provide a sharp view of the radio sky from from $\simeq 50$ MHz to $\simeq 14$ GHz – it is important to go beyond the conservative approach of Ref. [15], and to address in detail a number of important issues like the dependence of the bound on the mass distribution function, and the uncertainties inherent to the accretion physics.

The goal of this paper is precisely to go beyond the conservative estimates presented in Ref. [15]. We will first present a robust and updated upper bound on the DM fraction in PBHs that takes into account the possibility that PBHs have a general ‘extended’ mass function. We will then implement state-of-the-art numerical simulations of gas accretion onto moving BHs in order to more carefully characterize the expected multi-wavelength emission of PBHs, and to put the upper bound on the contribution of PBHs to the dark matter abundance itself on solid physical grounds.

2 Impact of the mass distribution on the astronomical bounds.

2.1 Introductory considerations

The expected mass distribution of PBHs has been a subject of much discussion since the first pioneering papers from Carr and Hawking (see, for instance, [16]). Although narrow mass distributions, well approximated by a Dirac delta, could be produced e.g. in the context of single-field inflation models [17, 18, 19, 20, 21], extended mass distributions (EMD) can arise from a variety of mechanisms: For instance, broad log-normal mass distributions peaked at large stellar masses can arise within hybrid inflation scenarios [22]; power-law distributions, on the other hand, can arise within inflation models featuring either a spectator field with a blue spectrum or a running spectral index [23]. A variety of other mechanisms, such as phase transitions [24, 25], or collapse of cosmic strings [26] or other topological defects, can also generate black holes with extended mass functions.

The impact of EMDs on constraints on the PBH abundance was first analyzed in [4, 27, 28]. An improved remapping method to convert the constraints was later proposed in [29] and applied to lensing and CMB constraints.

Given these considerations, it is very useful to analyze the impact of different types of EMDs on the radio and X-ray bounds reported in [15]. Following the approach of that paper, we consider a set of Monte Carlo simulations featuring a population of PBHs that trace the DM distribution in the Galaxy. The simulations include a realistic distribution of the interstellar gas, and rely on conservative assumptions regarding the physics of accretion and the subsequent non-thermal emission in the radio and X-ray bands: We refer to Ref. [15] for the details about the procedure and the astrophysical ingredients, here we recall the main features:

• $-$

We assume that the PBHs follow a Navarro-Frenk-White profile [30] (for a discussion of the impact of this assumption, see Sec. 4), and the gas distribution in the inner Galaxy is modeled as in [31].

• $-$

We assume that the PBH velocity distribution at a distance $R$ from the Galactic Center is a Maxwell-Boltzmann with mean velocity

[TABLE]

For the mass distribution in the Milky Way, $M(R)$, we consider the recent state-of-the-art model taken from [32].

• $-$

We conservatively assume, following, e.g., Ref. [33], (see Section 3 for a more complete discussion of the accretion physics) that the accretion rate of isolated PBHs is a small fraction $\lambda$ of the Bondi-Hoyle-Littleton rate,

[TABLE]

with,

[TABLE]

where $G$ is the gravitational constant, $M_{\mathrm{BH}}$ is the BH mass, $\rho$ is the density of the ambient gas, $v_{\mathrm{BH}}$ is the velocity of the BH (expressed here as a simple scalar quantity), and $c_{\mathrm{s}}$ is the sound speed of the gas. The adopted value $\lambda=0.02$ is consistent (albeit approaching the upper bounds) with isolated neutron star population estimates and studies of active galactic nuclei accretion [34, 35, 36].

• $-$

We model the radiative efficiency $\eta$ as follows:

[TABLE]

with $\eta=0.1\dot{M}/\dot{M}{\rm{}_{crit}}$ for $\dot{M}<\dot{M}\rm_{crit}$. The critical accretion rate $\dot{M}_{\mathrm{crit}}$ is defined in terms of the Eddington luminosity $L_{\mathrm{Edd}}$ as

[TABLE]

for which we fix the radiative efficiency to $\eta_{0}=0.1$, following, e.g., [33]. Since all the PBHs in our simulations fall below the critical accretion rate, they are inefficient accretors: As a consequence the luminosity scales non-linearly with the accretion rate: $L\propto\dot{M}^{2}$ (see, e.g., [37, 38] for more details about this regime).

• $-$

We compute the number of black holes in our simulations that would appear as point sources in the radio (in the GHz domain) and in the soft X-ray band (0.5-8 keV) above the detection threshold of the VLA and Chandra experiments respectively. To do so, we assume that $30\%$ the bolometric luminosity estimated by means of Eq. 2.4 is radiated in the soft X-ray band, and we convert the X-ray to radio luminosity by applying the empirical fundamental plane (FP) relation [39], which applies for a large class of compact objects in different mass ranges, from X-ray binaries all the way up to active Galactic nuclei. We then compare the predicted number of radio and X-ray sources associated with PBHs in our simulations to the number of observed ones, in a small region of interest centered on the Galactic center (see [15] for the full detail) dominated by a massive molecular structure known as the Central Molecular Zone (CMZ), and place the bound on the abundance of PBHs accordingly.

In the following section we discuss how the results in Ref. [15] change when adopting an extended PBHs mass distribution $f(M)$.

2.2 Log-normal mass function

We start by considering a log-normal distribution written as:

[TABLE]

where $M$ is the random variable, $\mu$ is the median mass, and $\sigma$ is the standard deviation of the distribution in the $log(M)$ space. The normalization is set to 1 so that the actual mass distribution is given by $dN/dM=N_{\rm tot}f(M)$.

The results associated with this distribution are visualized in Fig. 1. In principle, since the constraints depend on the number of detectable black holes (i.e., the black holes emitting radio or X-ray radiation above the detection threshold), and given that PBHs with mass $M>\mu$ ($M<\mu$) get higher (lower) fluxes with the implementation of an extended instead of a monochromatic (i.e., Dirac delta) mass distribution, the bound could get stronger or weaker depending on how the competing effects on the high- and low-mass tails of the distribution balance. In our case, given the fact that the constraints are driven by a tiny fraction of PBHs that—due to the low BH velocity and high ambient density—radiate above threshold, the implementation of log-normal mass distributions increases the fluxes of the PBHs close to the detection threshold, which in turn increases the number of detectable PBHs resulting in stronger constraints with increasing $\sigma$.

2.3 Power-law mass function

Let us now turn our attention to power-law EMDs. We parametrize the power-law EMD (normalized to 1) as in [29]:

[TABLE]

where the Heaviside step functions $\Theta$ implement the low- and high-mass cutoffs, and the normalization factor $\mathcal{N}$ is given by

[TABLE]

Such distributions have been considered since the early works by Carr and Hawking (see e.g. [40] and [26]), and is expected, e.g., if PBHs form out of the collapse of large density perturbations, or cosmic strings. We recall that the exponent $\gamma$ depends on the epoch of collapse. If the equation of state of the Universe at PBH formation is paramerized in the usual way as $P=w\rho$, the index $\gamma$ can be written as:

[TABLE]

For $w>-1/3$, as expected in an accelerating Universe, $\gamma$ covers the [-1, 1] range. Following [29], we have considered two reference values: $\gamma=-0.5$, which corresponds to PBH formation during the radiation-dominated era ($w=1/3$), and $\gamma=0$, which corresponds to PBH collapse during matter domination ($w=0$). The X-ray and radio bounds based on the power-law distribution for these relevant values of the power-law index $\gamma$ are visualized in Fig. 2. The color code in that plot corresponds to the PBH abundance excluded at $5\sigma$, for both the radio (left), and X-ray (right) case: the excluded value can be as low as $5\cdot 10^{-2}$ for large values of $M_{\rm max}$. We point out the very mild dependence on $M_{\rm min}$, and vice versa a strong dependence on $M_{\rm max}$. This behaviour is due to the fact that most of the constraining power actually comes from the high-mass tail, given the non-linear increase of the radiative output with $M$ and the subsequent significant increase of the number of high-mass PBHs emitting above threshold when $M_{\rm max}$ is increased.

2.4 A specific realization: PBH formation by vacuum bubbles.

Let us now consider in more detail, as a case study, a specific inflationary model that predicts a broad power-law spectrum of PBHs. Ref. [41] first noticed that non-perturbative quantum effects may play a relevant role during inflation. In particular, that paper proposed scenarios in which vacuum bubbles form during the inflationary phase due to the presence, in the underlying particle physics model, of vacua characterized by a lower energy density with respect to the one associated to the false vacuum that actually drives the inflationary expansion of the Universe. These bubbles are pulled inwards by the negative pressure of vacuum in its interior, and may collapse to PBHs after inflation. In ref [42], a broken power-law mass spectrum is derived by means of numerical simulations. This spectrum extends over many decades in mass, and—given the many existing constraints—the total density of PBHs in this scenario is bound to be below $10\%$ of the dark matter density; interestingly, the number of predicted PBHs in the range covered by Virgo and LIGO are compatible with the rate of binary-black-hole mergers reported by those collaborations.

Inspired by these results, we consider a mass distribution of the form:

[TABLE]

We also introduce a cutoff at $M_{\rm max}=1000$ $M_{\odot}$, given the very strong constraints on the PBH abundance in the high-mass domain (mainly from CMB measurements).The $5\sigma$ constraint is strong and very weakly dependent on the parameters involved: The maximum allowed DM fraction in the form of PBHs in such a scenario lies within the range $[5-6]\cdot 10^{-2}$.111For ease of comparison with [42], this range translates to $[4-5]\cdot 10^{-2}$ and $[7-9]\cdot 10^{-3}$ at $2\sigma$ significance for X-ray and radio, respectively, thus providing an improvement to the existing constraints by a factor of 2 to 15.

2.5 Remapping from monochromatic constraints.

It is now useful to consider the remapping method proposed in [29] and compare that prescription to our results, as a further validation of our computations. The remapping procedure is based on the principle that different PBHs contribute independently to the observable of interest. Therefore, the observable of interest is, in general, proportional to a weighted integral over the mass distribution $d\Phi/dM$,222To follow the notation of [29] and avoid confusion with the PBH dark matter fraction $f_{\mathrm{PBH}}$, here the mass distribution is denoted by $d\Phi/dM$ instead of $f(M)$. The meaning is exactly the same as before.

[TABLE]

with the function $g(M,\{p_{j}\})$ encoding the underlying physics that determines the contribution of PBHs of mass $M$ to the observable. In our case, the observable of interest is the number of accreting black holes radiating above the detection threshold, with $g$ containing all the details regarding the physics of accretion.

Since Eq. 2.11 holds for any mass function, one can compare to the case of a Dirac delta and write:

[TABLE]

where the superscript ${\rm MMD}$ refers to the monochromatic mass distribution. This equation implies that there exists a MMD, peaked at a given equivalent mass $M_{eq}$, that gives rise to the same value of the observable as the EMD under consideration: Therefore, the value of the constraint for an EMD can be read off from the delta-function constraint computed at $M_{eq}$.

We calculated the function $g(M,\{p_{j}\})$ for the X-ray and radio bound (see Appendix A for the details), and compared the converted constraints to our results obtained by means of Monte Carlo simulations. The comparison is shown in figure 3, the agreement is very good, and provides a solid, independent validation of our results.

3 Improved treatment of gas accretion onto PBHs: The impact on the bounds

The uncertainties associated with the physics of accretion play a major role in the estimation of the number of PBHs emitting above the threshold flux of successful detection. The constraints presented in [15] were obtained based on the principles of weak black hole accretion, derived from our current understanding of black holes observed to accrete at very low rates, e.g., Sagittarius A*, the black hole at the Galactic center (see, e.g., [43, 44] for reviews on this black hole). Such low-level accretors lead to models for radiatively inefficient accretion flows (RIAFs), e.g. [37], and the non-linear scalings of radiative efficiency discussed in Section 2.1 invoke these accretion scenarios. The accretion formalism of [15] makes use of this radiative inefficiency scaling within the framework of Bondi-Hoyle-Littleton accretion. In this setup the accretion rate onto the isolated PBHs was parameterised by $\lambda=\dot{M}/\dot{M}_{\mathrm{BHL}}$, and this incorporated all our ignorance regarding the energetics of the gas as it travels from the Bondi radius to the inner regions of the accretion flow. Ref. [15] then separately characterized the efficiency of converting the available gravitational energy into radiation through $\eta$, and factored in a two-phase gas regime, adjusting the sound speed according to whether or not the gas is likely ionized. Whilst this provides a reasonable framework within which to estimate the BH luminosity, a more physically-consistent approach would involve a more self-consistent estimate for how the accretion rate depends on the velocity of the black hole in the presence of radiative feedback.

In the following sections we present an improved treatment of accretion onto our population of PBHs with radiative feedback.

3.1 Radiation-hydrodynamic simulations

In this section, we place our bound on more solid ground by implementing the results of a set of numerical simulations, presented in Ref. [45], aimed at characterizing the Bondi–Hoyle-Lyttleton problem—the accretion of baryonic matter onto an isolated, moving compact object—in the presence of significant radiative feedback. These results build on previous findings discussed in Ref.s [46, 47], and are based on non-relativistic hydro-dynamic simulations that account for photo-heating and photo-ionization of the environment surrounding the compact object, due to the emission of UV and X-ray photons by the accretion disk and subsequent formation of an ionization front.

3.2 Analytic formalism

The fundamental distinction between the previous approach of simply parameterising the accretion rate in terms of some ratio of the Bondi-Hoyle-Littleton accretion rate, $\dot{M}_{\mathrm{BHL}}$, is that here we can instead prescribe the effective accretion rate within the ionized region behind the shock generated due to the motion of the BH. Thus here we may re-write the Bondi accretion formulae in terms of physical characteristics within this ionization region.

Taking first the basic definition of $\langle{\lambda_{\mathrm{B}}}\rangle=\langle{\dot{M}}\rangle/\dot{M}_{B}$ as the mean accretion rate with respect to the Bondi rate (thus in contrast with the definition used by [15] and shown in Equation 2.2, such that $\langle{\lambda_{\mathrm{B}}}\rangle=\lambda\dot{M}_{\mathrm{BHL}}/\dot{M}_{\mathrm{B}}$), Ref. [45] finds

[TABLE]

where the subscript ‘in’ denotes the relevant quantities within the ionized region, while $\rho$ and $c_{s}$ are the ambient medium density and sound speed, respectively, as in Eq. 2.3.

Now we switch to dimensionless notation such that the Mach numbers, temperatures, and densities are given by $\mathcal{M}_{\mathrm{in}}\equiv v_{\mathrm{in}}/c_{\mathrm{s,in}}$ and $\mathcal{M}\equiv v_{\mathrm{BH}}/c_{\mathrm{s}}$, $\Delta_{T}\equiv T_{\mathrm{in}}/T$ and $\Delta_{\rho}\equiv\rho_{\mathrm{in}}/\rho$, respectively. Since the gas within the ionized region is assumed to be fully ionized Hydrogen, the mean molecular weight of the gas decreases by a factor of 2, and thus $c^{2}_{\mathrm{s,in}}/c^{2}_{\mathrm{s}}=2\Delta_{T}$. The accretion rate is thus

[TABLE]

where we have adopted the Bondi formula for an isothermal spherically accreting stationary BH, with the factor $e^{3/2}$ coming directly from the assumption of an isothermal gas.

In order to calculate the accretion rate given by Equation 3.2 we need explicit expressions for two terms: 1) the Mach number within the ionized region, $\mathcal{M}_{\mathrm{in}}$, and 2) the density ratio across the ionization shock front $\Delta_{\rho}$.

The Mach number within the ionized region is related to the ambient Mach number by

[TABLE]

Regarding the latter term, we recall that—as detailed in [45]—the radiation-hydrodynamic simulation results are interpreted through the transitions between two different types of ionization fronts as the BH mach number increases: rarefied, or $R$-type, and dense, or $D$-type. The key difference lies in the contrast in gas density across the I-front: an $R$-type I-front shows rough equilibrium between these two densities, $\rho_{\mathrm{in}}\sim\rho$, whereas a $D$-type front exhibits a lower density downstream than the upstream density, $\rho_{\mathrm{in}}\leq\rho$. The full expressions as a function of the ambient Mach number and temperature ratio, which we will not derive here, are given by

[TABLE]

Here the Mach number boundaries for $D$-type and $R$-type I-fronts are given by $\mathcal{M}_{D}=\sqrt{2\Delta_{T}}-\sqrt{2\Delta_{T}-1}$ and $\mathcal{M}_{R}=\sqrt{2\Delta_{T}}+\sqrt{2\Delta_{T}-1}$ respectively. Substituting Equations 3.4 and 3.3 into Equation 3.2 allows us to fully characterize the PBH accretion rates in accordance with the numerical simulation results of [45]. Motivated by empirical values of the temperature in the CMZ (see [48] and references therein, in particular [49]), we adopt the same temperature assumptions as in [15] for the neutral and ionized regions respectively, $T=10^{2}$ K ($c_{\mathrm{s}}=1$ km/s) and $T_{\mathrm{in}}=5\cdot 10^{3}$ K ($c_{\mathrm{s,in}}=10$ km/s), giving $\Delta_{T}=50$.

The accretion rate as a function of BH velocity is shown in Figure 4 (solid blue line), where we have compared the rates to the previous basic accretion assumptions of [15] (dashed and dash-dotted lines). The plot clearly shows that the key feature of the new approach presented above and based on [45] is the behaviour of the accretion rate as a function of the Mach number of the isolated moving black hole, captured by Eq. 3.2.

Let us comment further to convey a clearer physical intuition of this trend. While the classical Bondi–Hoyle–Lyttleton formula adopted in the previous work implies a simple, monotonic decrease of the accretion rate with increasing BH velocity, the simulations show a more complicated phenomenology. If the BH velocity is supersonic, but the Mach number is below a critical value $\mathcal{M}_{R}$ ($\mathcal{M}_{R}\simeq 4$ for $T\simeq 10^{4}$ K), a dense bow shock forms in the upstream region; behind the bow shock, a D-type (dense) ionization front develops, and a cometary-shaped HII region can be identified, characterized by low density and velocity. In this regime, the gas velocity in the reference frame of the moving BH decreases with increasing BH velocity, hence the accretion rate follows the opposite trend with respect to the Bondi-Hoyle-Lyttleton formula, i.e. it increases with increasing BH velocity (see Fig. 4). Conversely, when the BH velocity is above $\mathcal{M}_{R}$, the ionization front becomes R-type (rarefied), and the accretion rate decreases with the BH velocity, effectively returning the accretion process to Bondi-Hoyle-Littleton-like accretion.

3.3 Impact on the PBH constraints

The accretion formalism described above has a noticeable impact on the constraints derived in [15]. Upon further inspection of Fig. 4, it is possible to build an intuitive understanding for the constraints as a function of PBH mass and velocity. With the new prescription (blue line), the accretion rates of PBHs with low velocities are significantly suppressed and should no longer contribute to the constraint; instead, the accretion rate of PBHs with higher velocities is enhanced and therefore these objects are more likely to radiate above threshold.

Our simulations actually confirm this intuition and clearly show that the bound no longer originates from the PBHs at the low velocity tail, but from the PBHs with velocities around $2c_{s,\mathrm{in}}=20$ km/s, i.e. near the peak of the accretion rate as plotted in Fig. 4.

The resulting constraints are shown in Fig. 5: They are two orders of magnitude stronger than those presented in [15], excluding PBH DM fractions above $\sim 10^{-3}$. This remarkable effect is mainly due to our result being sensitive to higher velocities—as pointed out above—which correspond to a larger fraction of the PBH population.

As far as the mass dependence is concerned, we notice that the radio and X-ray bound have different behaviors. This is due to two competing effects: 1) On the one hand, the accretion rate rises with increases mass, increasing the probability of detection; 2) on the other hand, the number density of PBHs decreases with increasing mass given the fixed reference dark matter mass density.

For the radio bound, the former effect dominates since the constraining power is actually driven by a very small population of PBHs radiating close to the detection threshold, while the number of observed BH candidates in the radio is 0 as discussed in [15]: therefore, the radio bound is not present for masses lower than 20 $M_{\odot}$, while for larger masses it “kicks in” and gets stronger with increasing mass.

Concerning the X-ray bound—given the higher sensitivity of Chandra and the resulting large population of PBHs radiating above threshold—the latter effect dominates: hence, the X-ray bound gets progressively weaker with increasing PBH mass.

Finally, let us point out that the estimated radio luminosities of the PBH population, relative to their X-ray luminosity, also increases with BH mass in accordance with the FP scaling relation [39], which also contributes to the relative strengthening of the radio bound at higher BH mass with respect to the X-ray bound.

4 Discussion

As a first discussion point, we want to emphasize the role of several critical aspects of accretion physics. As pointed out in [47] and [45], the flow of gas being accreted onto a moving, isolated compact object features a rich and complex phenomenology.

An important feature is time dependence: According to simulations, at intermediate densities ($n_{H}\sim 10^{3}$ cm*-3*, which is a close to the values measured in the CMZ clouds) the accretion rate shows an oscillatory behavior, and intermittent bursts can occur. These bursts are caused by an instability in the ionization shell in the upstream direction, and are seen to occur in radiation-hydrodynamic simulations in the Mach number range $2.5\leq\mathcal{M}\leq\mathcal{M}_{R}$. The instability causes the shell to fragment and reform intermittently, with explosive accretion occurring while the shell is fragmented, and the reformation occurring in roughly the sound crossing time of the HII region. These bursts can increase the signal by a factor of 3 or more for a relatively short time (on the order of 1000 years, with a recurrence time of   5000 years, both timescales depending on the gas temperature and density), and can significantly increase the probability of detecting a moving PBH that is crossing a dense cloud: we chose to neglect this effect, given the conservative approach we have followed, and the difficulty and the uncertainties involved in modeling this effect; however, it is useful to consider the fact that the bound can get even stronger if this behavior is taken into account.

In addition, we recall that, while we are considering isolated black holes accreting gas as they traverse the CMZ, most of our knowledge on accreting BHs comes from those that are situated in X-ray binary systems: these are the closest such systems we can actually observe and study in detail. In this work, we have neglected the possibility that the more well-known outbursting behaviour of X-ray binaries could occur due to a critical mass forming in a disk-like structure during accretion (i.e. the disk instability model, see, e.g., [50]). As mentioned by [45], the formation of an accretion disk-like structure akin to those found in the many accreting sources is not resolved in their simulations. As such we cannot be certain that such structures will not form around accreting isolated PBHs, and that this would not lead to the thermal instabilities expected in dense, thin accretion disks. The possibility of isolated BHs accreting from the ISM exhibiting such transient behaviour has been considered already [51, 52]. If density perturbations within the CMZ exist on similar size scales to the Bondi sphere, the specific angular momentum of the gas may be sufficient for the gas to circularize prior to accreting [53]. Whilst [52] suggest that gas densities in the accretion flows of isolated BHs can become sufficient to trigger Hydrogen ionization instabilities (and thus transient outbursts), as shown by [51], the lack of detection of transient isolated systems (be they PBHs, astrophysical BHs, or Neutron stars) makes it likely that only a small fraction of isolated BHs actually exhibit outbursts of this nature. We have neglected this possibility, but note that folding in the predicted outburst rates of PBHs in the CMZ may have the effect of increasing the strength of the PBH DM bounds.

Another key assumption behind the constraints derived here and in [15] is the presence of jets launched from the accretion flows of the PBHs. As discussed by, e.g., [54], radio searches for isolated BHs should be more successful than X-ray searches, since the jet kinetic power dominates the total accretion power with decreasing luminosity. This is a direct result of the observationally determined radio and X-ray flux correlations of X-ray binaries ($L_{R}\propto L_{X}^{0.7}$; [55, 56, 57, 58]), and is shown to be a mass-scaled relation [59, 38, 60], ultimately leading to the FP relation we use here to estimate the radio luminosities of our population of accreting PBHs [39]. As shown in Fig. 5, whilst the X-ray bounds on the PBH DM fraction are the most constraining at lower masses, at higher BH masses the radio bounds become stronger. This is due primarily to the vast difference in detected sources in the respective Chandra and VLA Galactic center surveys, but is is also due in part to the mass-scaling and the dominance of jet kinetic power at the lowest accretion rates. With future radio telescopes such as SKA, we should expect the bounds to be dependent mostly on the radio surveys. Here, just as was done by [15], the predicted radio fluxes are based on the FP relation [39]. In order to apply the FP relation, we have to assume the presence of self-absorbed jets in our systems, which are emitting GHz radio waves with an almost flat, optically thick spectrum. We remark that only the radio bound depends on this assumption, since the radio flux is estimated by means of a conversion from the X-ray band based on the FP relation. Ref. [45] discuss the fact that the dynamical and radiative implications of the presence of jets is not considered in their simulations.

The prevalence of jets in accreting BH systems across the mass scale at the lowest accretion rates (see, e.g., [39, 58]) suggests that the structure of accretion flows (optically thin and geometrically thick) associated with low accretion rates are ideal for the launching of jets. This has been explored in the case of isolated BHs [61, 54], whereby the assumed conditions necessary for jet launching include a spinning BH and the presence of a magnetic field. We do not explore the likely distribution of PBH spins here, nor predict the strength of the magnetic field in the accretion flow, but we note that such conditions would be necessary. Thus the question of how this would impact the radiation-hydrodynamic simulations of [45] is a caveat to our implementation in the simulations, but only for the derived radio bounds.

Another relevant source of uncertainty is the phase-space distribution of both the visible and the dark matter in the inner part of the Galaxy.

Regarding the former, as mentioned e.g. in [62, 48] and references therein, an uncertainty of a factor of $\simeq 2$ or more about the total mass of the CMZ has to be taken into account.

As far as the latter is concerned, although there is evidence for the presence of dark matter inside the solar circle [63], the actual distribution at small Galactocentric radii is very uncertain. We assumed the NFW profile as a benchmark in this work, but we tested how the bound changes using different profiles. In particular, we remark that the scenario of PBHs composing all dark matter (i.e. $f_{\mathrm{DM}}=1$) is still ruled out at the $5\sigma$ confidence level even under the extreme assumption that the dark matter profile has a core as large as the size of the Galactic bulge.

As a further validation of our results, we considered an improved estimate of the velocity distribution based on the Eddington formalism [64]. Under the assumption of spherical symmetry and isotropy, the Eddington inversion formula allows us to compute the velocity distribution at different Galacticentric radii $f(r,v)$ given a density profile $\rho(r)$. This formalism has been successfully applied in the context of the determination of the DM density in the Solar neighborhood [65]. In our case, since the constraining power mainly comes from the very inner regions of the Galaxy, the assumptions behind the formula are put on even safer grounds. Moreover, while in the previous paper the constraining power in the momentum space mainly came from the very-low-velocity tail of the distribution, now—given the increase of the accretion rate with the PBH velocity—PBHs above the detection threshold cover a wide range of velocities. Although the impact of the Eddington formalism, which departs from the Maxwell-Boltzmann approximation especially at large speed, may thus be in principle significant, we have verified that it has little impact on the X-ray and radio bound.

Regarding the velocity distribution, we stress that the velocity we are interested in is the relative one between the black hole and the gas, which can be affected by the turbulent motion of the gas itself. This point is stressed by the authors of [66] who revisited the bound proposed in our previous work by adding the turbulent speed to the denominator of the accretion rate (i.e., to the denominator of Eq. 2.3, substituted into Eq. 2.2) as follows: $\dot{M}=4\pi\lambda(GM_{BH})^{2}\rho\left(v_{BH}^{2}+v_{turb}^{2}\right)^{-3/2}$. However, they have overestimated the impact of turbulent motion by adopting as a reference value for $v_{turb}$—taken from [48] and based on spectral line observations of molecules such as HNCO, N2H+, and HNC—an upper limit, not an estimate, of the speed of turbulent motion [48]. Even assuming the maximum value for $v_{turb}$, our results do not change significantly, since the bulk of the high-luminosity radio and X-ray sources, as we stressed several times above, do not arise from PBHs in the very-low velocity tail.

As a further caveat, let us mention that there has been a debate in the literature about the clustering level of PBHs [67, 68, 69, 70] and its impact on the bounds. Whilst a quantitative estimate of this effect is beyond the scope of the present work, we remark that a significant clustering could in principle go in the direction of strenghtening the bound, since a cluster of PBHs crossing a dense gas cloud would be more easily detected, and leave a significant and observable imprint on the environment due to the localized, high flux of ionizing photons from the accreted gas.

As a final point, we stress that PBHs searches in the radio are particularly promising, in view of the upcoming radio facilities, including the SKA observatory [71] and ngVLA [72], which are expected to provide a very significant increase in sensitivity. One of the main challenges will be then to devise effective multi-wavelength strategies to disentangle accreting PBHs from other radio sources, and in particular from the guaranteed population of astrophysical black holes, since the current bound is approaching the estimated density associated with the population of stellar-mass astrophysical black holes in our Galaxy that form through conventional stellar evolution. Such a population could amount to $\sim 10^{8}$, see e.g. [73, 33]: this estimate stems from a modeling the star formation history of the Milky Way and our knowledge about the local mass density of stellar remnants. Very little is currently known about this population, with a few exceptions (namely, the X-ray binaries—i.e. the black holes that are accreting mass from a companion star—and possibly a black hole candidate part of a wide binary system recently announced in [74]).

5 Conclusions

In this paper we have discussed astronomical searches for primordial black holes based on the analysis of radio and X-ray data from the Galactic center region. We have updated the contraints on the abundance of PBHs presented in [15] by taking into account the impact of an extended mass function and by introducing a more detailed treatment of the physics of gas accretion onto moving PBHs. We have shown that, in general, extended mass functions lead to constraints that are more stringent by a factor of a few with respect to the case of a monochromatic mass function. We have validated our results via a semi-analytical remapping procedure.

As for the physics of gas accretion, we have implemented, for the first time in this context, the results of a set of state-of-the-art numerical simulations that capture the rich phenomenology of gas accretion onto a moving compact object in the presence of radiative feedback. This new approach allowed us to obtain a strong bound on the PBH abundance, which is at the level of $\sim 10^{-3}$ of the dark matter density (and close to the expected level for the astrophysical black hole density). The bound comes from both X-ray and radio data, and with respect to the previous upper limit reported in [15], it appears less sensitive to the details of the low-velocity tail of the velocity distribution.

In conclusion, astronomical data strongly disfavor primordial black holes as the main constituent of the dark matter in the mass range under consideration, and can be considered a promising path of detection for both a sub-dominant population of those objects (with a number density within the allowed range), and the expected population of astrophysical black holes, either isolated or as part of wide binary systems, the detection of which is a current challenge of modern astrophysics.

Acknowledgements

We thank G. Ballesteros, J. Garcia-Bellido, and T. Maccarone for important comments and suggestions. D.G. has received financial support during the last stage of this project through the Postdoctoral Junior Leader Fellowship Programme from la Caixa Banking Foundation. R.M.T.C. acknowledges support from NASA Grant No. 80NSSC177K0515.

Appendix A Derivation of $g(M,\{p_{j}\})$

As mentioned in Sec. 2.5, the observable of interest is the number of PBHs radiating above the detection threshold. With each PBH contributing independently to this number, the observable is proportional to Eq. 2.11. Therefore, the number of sources above the detection threshold $N_{S}$, expressed as an integral over the flux distribution of the PBH population $dN/d\phi$, can be equated to an integral over the mass distribution,

[TABLE]

Here the proportionality between the observable and Eq. 2.11 has been fixed by the usage of the total mass distribution $dN/dM\equiv N_{\mathrm{tot}}\,d\Phi/dM$, with the total number of PBHs $N_{\mathrm{tot}}$ having an implicit dependence on the PBH dark matter fraction $f_{\mathrm{PBH}}$. Now, an intuitive description of $g$ becomes visible; it is the fraction of PBHs with mass $M$ that contribute to the number of sources with emission above the detection threshold.

In principle, the flux detected at Earth originating from an accreting PBH depends on many different parameters. Within the simulations, this is limited to three randomly chosen numbers: the PBH mass, the PBH velocity and its position. To obtain the flux distribution, one needs to convert the mass, velocity and position distributions following the calculations performed in the simulations. Given the complexity introduced by the random locations, we proceed only with a random velocity and mass. Therefore, we are making the assumptions of a constant velocity distribution, distance and gas density.

For simplicity, we also make the assumption of a bijective mapping of the BH velocity and mass to its radiative flux. Clearly, this assumption is only applicable to the accretion rate considered in Sec. 2, and not to that considered in Sec. 3 where multiple velocities can correspond to the same accretion rate (see Fig. 4). One could in principle extend the following formalism also for Sec. 3 by introducing thoroughly chosen breaks in the velocity integration limits, but this is beyond the scope of this paper. Thus the remainder of the Appendix only applies to Sec. 2.

The aforementioned assumption allows us to perform a change of variables on the flux distribution in Eq. A.1. By consequently comparing it with the right-hand side, an implicit expression for $g$ can be obtained,

[TABLE]

Here $p$ is the normalized velocity distribution, and the function $h(\phi,M)$ maps a given flux $\phi$ and mass $M$ to the velocity required for a PBH with mass $M$ to obtain this flux. Given that under our assumption the flux is a monotonically decreasing function of the velocity, $g$ is nothing more than the CDF of the velocity distribution, properly taking into account the detection threshold. Since we are considering a Maxwell-Boltzmann velocity distribution, a more explicit version of $g$ is given by its CDF,

[TABLE]

with the mass dependence hidden in $h_{\mathrm{thres}}\equiv h(\phi_{\mathrm{thres}},M)$. The minimum mass $M_{\mathrm{min}}$ is introduced due to the fact that in the case of a steady black hole the accretion is limited by the sound speed of the accreting gas, and thus light black holes may never reach fluxes above the detection threshold.

The last step is to find an explicit expression for the function $h$, containing the calculation of the flux from the mass and velocity as in the simulations. This depends on the flux of interest and are given by

[TABLE]

for the 1.4 GHz radio and $0.5-8$ keV X-ray flux respectively. Here $c_{s}$ is the sound speed of the surrounding medium, and $X=0.0448\,(\mathrm{mJy})^{1.45/3}\,(M_{\odot})^{-3.88/3}\,(\mathrm{km/s})^{2}$ and $Y=0.0277\,(\mathrm{ph}\,\mathrm{cm}^{-2}\mathrm{s}^{-1})^{1/3}\,(M_{\odot})^{-1}\,(\mathrm{km/s})^{2}$ are constants.