Multi-wavelength detectability of isolated black holes in the Milky Way

Isolated black holes in our Galaxy have eluded detection so far. We present here a comprehensive study on the detectability of isolated stellar-mass astrophysical black holes that accrete interstellar gas from molecular clouds in both the local region and the Central Molecular Zone. We adopt a state-of-the-art model for the accretion physics backed up by numerical simulations, and study the number of observable sources in both the radio and X-ray band, as a function of a variety of parameters. We discuss in particular the impact of the astrophysical uncertainties on our prediction for the number of bright X-ray sources in the central region of the Galaxy. We finally consider future developments in the radio domain, and assess the potential of SKA to detect a population of astrophysical black holes accreting gas in our Galaxy.


INTRODUCTION
A large population of Stellar-mass black holes (BHs) is believed to exist in our Galaxy: several works (Shapiro & Teukolsky 1983;Samland 1998;Caputo et al. 2017) agree on an order-of-magnitude estimate of ∼ 10 8 BHs produced by the collapse of massive stars in the Galaxy. This large number of BHs has eluded detection so far, with the exception of 60 objects 1 that are part of binary systems and accrete a significant amount of mass from a companion star: such systems are known as black-hole X-ray binaries (BH-XRBs) and appear as bright X-ray sources in the sky.
Here, we focus on isolated BHs, whose detection represents a major challenge in modern astronomy. The problem has been discussed in several recent studies focused on the future detectability of the multi-wavelength radiation possibly emitted by these objects when they accrete matter from interstellar clouds either in the vicinity of the Sun (Maccarone 2005;Fender et al. 2013) or in the regions near the Galactic centre, as recently discussed in Tsuna & Kawanaka 2019;Tsuna et al. 2018. The latter region seems particularly promising. In fact, an extensively studied complex of giant molecular clouds usually called Central Molecular Zone (CMZ) ★ E-mail: francesca.scarcella@uam.es † E-mail: daniele.gaggero@uam.es 1 http://www.astro.puc.cl/BlackCAT/ characterize a vast region that extends up to 200 pc away from the Galactic Centre, and provides a huge reservoir of molecular hydrogen. The high-density clumps within these region appear as ideal targets for astronomical searches of compact objects based on gas accretion. In order to perform these analyses, it is crucial to develop an accurate modeling of (i) the accretion process of baryonic matter onto a compact object, in order to compute the expected accretion rate as function of the environmental conditions and the BH mass and speed; and (ii) the non-thermal spectrum of the radiation emitted by the accreted gas in different wavelengths, from radio all the way up to the hard X-ray band.
Regarding the former aspect, a simplified model for the accretion physics is adopted in the aforementioned papers, based on a rescaling of the Bondi-Hoyle-Lyttleton (BHL) formula. However, the accretion problem exhibits a richer phenomenology. Radiation feedback plays a crucial role, and is expected to shape the surrounding environment in a remarkable way, forming a cometary-shaped ionized region around the compact object. A state-of-the-art treatment of accretion in the presence of radiative feedback, backed up by numerical simulations, was presented in a series of works (Park & Ricotti 2011, 2012Sugimura &Ricotti 2013, PR13 hereafter) and showed a more complicated behaviour of the accretion rate with respect to the BH speed compared to the BHL scaling. This model features an increase at low velocities, and an asymptotic decrease to the BHL rate at large velocities. The PR13 results were recently applied to compact object searches in the context of the quest for BHs of primordial origin (Manshanden et al. 2019). However, a comprehensive study of this model related to the population of astrophysical BHs is still lacking.
As for the prescription for the multi-wavelength (radio/X-ray) non-thermal emission, the absence of any observational detection of isolated accreting stellar-mass BHs means any physical model one can develop carries large uncertainties. However, one can utilize the vast knowledge gathered in the studies of both accreting BH-XRBs and accreting supermassive BHs to inform such a prescription. In particular, the Fundamental Plane relation between radio and Xray emission (Plotkin et al. 2012a) derived from the study of the aforementioned systems can be a precious tool to characterize the relation between radio and X-ray emission.
The aim of this paper is to characterize the detection prospect of a population of astrophysical black holes in view of the future development of radio and X-ray astronomy (with particular reference to the SKA experiment for the former), taking into account the state-of-the-art accretion formalism from PR13.
The paper is organized as follows: In Secs. 2 and 3 we review the main aspects of accretion physics with particular focus on the PR13 model, and discuss the emission mechanisms that are most relevant in the radio and X-ray domain. Then, we present our methodology in Sec. 4 and 5 and turn our attention to the astrophysical black holes in the vicinity of the Sun in Sec. 6 estimating the number of detectable X-ray sources. We then focus on the Central Molecular Zone in Sec. 7, and present a comprehensive parametric study of the expected X-ray emission from astrophysical BHs in this region. Finally, in Sec. 8, we present a forecast about the number of radio sources potentially detectable by SKA in association with the population of isolated astrophysical BHs accreting gas in the CMZ.

SETTING THE STAGE: STATE-OF-THE-ART ACCRETION PHYSICS
In this section we describe in detail the accretion model introduced in PR13 and compare it to the textbook Bondi-Hoyle-Lyttleton model.

The Bondi-Hoyle-Lyttleton model
According to the BHL accretion model (Hoyle & Lyttleton 1939;Bondi & Hoyle 1944), the rate of accretion onto an isolated compact object of mass moving at a constant speed v BH is given by: where and s are, respectively, the density and the sound speed that characterize the ambient medium and is the gravitational constant. This formula is usually re-scaled by introducing a suppression factor . Fender et al. (2013) for instance, found that values larger than = 0.01 are excluded, under realistic assumptions, by observations of the local region, where a significant population of isolated BHs should be present. Tsuna & Kawanaka (2019) and Tsuna et al. (2018) consider values between = 0.1 and = 10 −3 . The parameter may effectively capture the outflow of material that is expelled from the Bondi sphere, as suggested by several authors (Blandford & Begelman 1999). Its introduction is also supported by the non observation of a large population isolated neutron stars in the local region (Perna et al. 2003) and the studies of nearby AGNs (Pellegrini 2005) as well as the supermassive BH at the Milky Way centre, Sagittarius A* (Baganoff et al. 2003). In conclusion, it seems that the BHL accretion formula may overestimate the accretion rate by orders of magnitude, although the physical mechanism behind this deviation is still disputed. One such mechanism is introduced in PR13, as we discuss below.

The Park-Ricotti model
The hydrodynamical simulations performed in PR13 show that, when radiative feedback is taken into account, a cometary-shaped ionized region is formed around the BH as it moves through the interstellar medium.
The model proposed in the same work combines the BHL formula with the modelling of the ionization front to obtain an analytical formula in agreement with the results of the simulations. We review it here. First, BHL accretion is assumed to hold within the ionized region, so that the PR13 accretion rate can be written as: where s,in and in are the sound speed and density of the ionized medium, and v in is its velocity relative to the BH. The sound speed s,in is a free parameter of the model (of the order of few tens of km/s). Now Euler's equations can be applied to the ionization front to express in and v in in terms of the corresponding quantities referred to the external neutral medium: its sound speed , its density and the relative velocity of the BH, v BH .
We briefly discuss here the derivation of these relations. Assuming a one-dimensional steady flux, the jump conditions obtained by applying Euler's equations for mass and momentum conservation at the ionization front are: where to obtain the second equation we have further assumed an ideal gas. These equations have the following solutions: In the high-velocity range, v BH ≥ v R , the correct solution is given by − in , while for v BH ≤ v D , the front is best described by by + in (Park & Ricotti 2013). In the intermediate-velocity range, v < v BH < v , no real common solution can be obtained for Eqs. 3 and 4. Correspondingly, as velocity lowers below v , the pressure wave building behind the ionization front detaches from it as a shock into the neutral material. Part of the flux of matter is deviated and the velocity of the gas beyond the shock is reduced to below v . This implies on the one hand that Eq. 3, stating the mass conservation through the front along the direction of displacement of the BH, is no longer valid. On the other hand, since the velocity past the shock is now lower than v , the jump conditions at the ionization front can be solved. One should now consider, however, two sets of jump conditions associated to the two fronts: the shock and the ionization front. Park and Ricotti instead observed through simulations that, in this regime, the velocity inside the ionized region is v in ≈ s,in . This relation, promoted to an equality, can be used together with Eq. 4 to compute the density in . This way, for v ≤ v BH ≤ v , we obtain: Thus we have in summary: and Plugging these equations back into Eq. 2 finally gives the desired accretion rate expressed in terms of the BH speed and the properties of the neutral medium it is moving through.
Notice in particular that in the velocity range v D ≤ v BH ≤ v R we get: which increases quadratically with the BH velocity, in sharp contrast to the behaviour of the BHL rate, which decreases with velocity. This behaviour is the main feature introduced by the PR13 model, and is due to the formation around these objects of the aforementioned bow shock that deflects part of the gas away from the BH. The velocity dependence of the BHL rate is recovered in the high velocity regime, v BH > v R , where both rates present a ∝ v −3 BH dependence. The complete velocity dependence of the PR13 rate is shown in figure 1, for varied gas densities, BH masses, and sound speeds of the ionized region. For comparison, the BHL rate with = 1 is also shown. We can observe in this figure how the BHL rate decreases monotonically with velocity, whereas the PR13 rate peaks at v BH = v R = 2 s,in and is suppressed at lower velocity by the presence of the bow shock.
For < v D , the rate increases again. However, notice that this transition typically happens at velocities of ≈ 0.1 km/s (see Eq. 6), which are of little relevance for this work and not shown in figure.
The different velocity dependence of the PR13 rate compared to the BHL rate has important consequences when studying the emission properties of a BH population characterized by a given velocity distribution. According to the BHL prescription, the lowvelocity tail of the population is the easiest to detect. Following PR13, the highest emissions are instead associated to BHs with intermediate velocities.
Furthermore, BHL can predict very high accretion rates if  Figure 1. Accretion rate as a function of the BH speed. We show the accretion rate obtained according to the PR13 and BHL models, as a function of the BH speed and other relevant parameters. For the BHL rate, we set the suppression factor = 1, to allow for a more direct comparison. The to models agree at high velocity, but predictions differ by many orders of magnitude in the low velocity range. the speeds are low enough, which can easily lead to overshooting experimental bounds, while the highest rates predicted by PR13 are orders of magnitude smaller.
As a final remark, using Eq. 10 to express the peak of the PR13 rate in terms of the Eddington accretion rate Edd : one can see that the accretion rate will always be highly sub-Eddington in the range of BH masses BH , gas densities and ionized sound speeds ,in we consider in this work. This is taken into account when defining the prescription for the luminosity in the next section.

EMISSION MECHANISMS
In order to translate the predicted accretion rates of a population of isolated BHs into a prediction of detectable sources we need estimates of the associated bolometric luminosity of a given source, and the Spectral Energy Distribution (SED). Our focus is on the Xray and radio luminosity, X and R , of the accreting BH. We can estimate both by considering the SED of the BH at varying accretion rates. We are particularly interested in emission from highly sub-Eddington accreting BHs (as predicted by Eq. 11). Therefore in this section we describe our simple framework of the properties of sub-Eddington flows onto isolated BHs based on our prior understanding of known weak accretors in nature: Galactic BH-XRBs, and lowluminosity AGN.
Radiatively Inefficient Accretion Flow (RIAF) models have spearheaded many studies regarding the emission mechanisms associated with such sub-Eddington accretion flows, with a focus both on Galactic BH-XRBs at low accretion rates, and the supermassive BH in the Galactic center, Sgr A* (Narayan & Yi 1994;Yuan et al. 2003).
The regularly cited RIAF models which first considered the emission processes in inefficiently radiating accreting BHs are classed as advection-dominated accretion flows (ADAFs; Narayan & Yi 1994;Esin et al. 1997). In the ADAF model, Bremsstrahlung radiation dominates the observed spectrum at the lowest accretion rates, and thus the spectrum has some curvature and its peak energy depends on the thermal gas temperature. As the accretion rate and gas density increase, Inverse Compton (IC) scattering begins to dominate (likely via SSC, i.e., Synchrotron photons become scattered), though if densities remain low enough, the spectrum will still show come curvature. At higher accretion rates, the IC scattering process is more efficient, resulting in a power-law-like spectrum in the X-ray band.
The key thermodynamic property of RIAFs is the inefficient cooling of ions due to the Coulomb decoupling of electrons and ions in the accreting, low-density plasma. Due to this decoupling, such flows are likely well described by a two-temperature electronion plasma, with only electrons radiating via Bremsstarhlung, Synchrotron, and inverse Compton scattering Esin et al. (1997). The models, as well as a plethora of observational evidence, show that such flows display radiative inefficiency, with = 0.1 / crit , where is defined by: and crit is the accretion rate below which we have a RIAF. We assume crit is the accretion rate corresponding to 1% of the Eddington luminosity with Edd = 0.1 (Fender et al. 2013).The bolometric luminosity in such inefficient states is thus given by thus exhibiting a quadratic scaling ∝ 2 (Esin et al. 1997), as opposed to the linear behaviour observed at higher efficiency.
A significant fraction of this bolometric luminosity is typically assumed to fall in the X-ray band. The / fraction is usually assumed to be 30% (Fender et al. 2013).
Radiatively inefficient Galactic BH-XRBs also almost ubiqitously launch outflows. Such outflows are typically in the form of steady, self-absorbed relativistic jets (Fender 2001). Ballistic, transient jets are also observed during spectral state transitions at higher accretion rates. Such jets are typically identified via their radio emission, and are shown to be present during quiescence, at luminosities below 10 −8 Edd (Plotkin et al. 2015). These steady jets exhibit a flat-to-inverted spectrum from radio through to IR frequenciesa consequence of self-absorbed synchrotron emission through the optically-thick regions of the jet (see, e.g., Markoff et al. 2005).
In addition, multiple studies have now established that the luminosity of such radio-emitting jets scales in a consistent and predictable way with the X-ray emitting plasma in the inner regions of the flow. This scaling, R ∝ 0.7 x (Corbel et al. 2000(Corbel et al. , 2003Gallo et al. 2003;Corbel et al. 2008;Miller-Jones et al. 2011;Gallo et al. 2014), is known as the radio-X-ray correlation, and applies to Galactic BH-XRBs (though other compact objects have been tracked in this phase space, showing similar but distinct trends of their own). Therefore, by using this analogy between isolated BHs and their low-accretion rate counterparts in binaries, we can infer that isolated BHs will similarly launch these steady jets, and thus their radio and X-ray luminosities will scale according to this radiatively inefficient track we observe in Galactic BH-XRBs.
However, in order to invoke a mass scaling and capture lowluminosity accretion onto a mass-variable population of isolated BHs, one must refer to the established connection between Galactic BH-XRBs and AGN-the Fundamental Plane of Black Hole Activity (FP; Merloni et al. 2003;Falcke et al. 2004;Plotkin et al. 2012b). The FP is, in a sense, an extension of the radio-X-ray correlation discovered in Galactic BH-XRBs to their supermassive counterparts. The FP is an empirical, parameterized relation between the X-ray luminosities, radio luminosities, and masses of hard state Galactic BH-XRBs with steady jets, and AGN of types which display similar X-ray emission characteristics. These include low-luminosity AGN (LLAGN), low-ionization nuclear emisssionline regions (LINERS), Faranoff-Riley type I (FRI) and BL Lacs. It is understood that the fundamental connection between these types of AGN and Galactic BH-XRBs is simply their sub-Eddington accretion rates, and the presence of radio-emitting jets. All sources in the FP display something resembling a power law spectrum in the X-ray band, and a flat-to-inverted radio spectrum. Thus, by invoking this scaling relation, and assuming that an isolated population of low-luminosity accreting BHs will display similar spectral properties, one can adopt the FP relation to scale their X-ray and radio luminosities with BH mass. We therefore adopt the following scaling relation, determined in one of the most recent such FP studies (Plotkin et al. 2012b): Where is the X-ray luminosity in the 2-10 keV band and is the radio luminosity at 5 GHz. Thus we have a simple prescription for both the X-ray luminosity, X , and radio luminosity, R , of an isolated accreting BH.

ESTIMATING THE NUMBER OF VISIBLE SOURCES
Having described our prescription for the accretion rate and the associated non-thermal emission, we now turn at assessing the number of isolated astrophysical black holes that are potentially detectable by the current (and forthcoming) generation of X-ray (with particular focus on the hard X-ray band) and radio experiments. In this section we summarize the main points of our methodology. Our main observable is the number of sources associated to a radiation flux above detection threshold, i.e. that satisfy > * , where is the flux at Earth defined as: = /(4 2 ) with r the distance of the source from Earth, and * is the threshold value.
Previous studies (Fender et al. 2013;Tsuna et al. 2018; Tsuna & Kawanaka 2019) obtained their estimates through Montecarlo simulations of a BH population. Here we adopt instead a semi-analytical approach which allows us to perform a comprehensive parametric study associated to the physical model of accretion physics described above.
Here we describe the general prescription for computing the number of visible sources, which we will apply with some variation to the two region of interest for this work: the local region (section 6) and the central molecular zone (CMZ, section 7). We obtain the probability for a BH to emit above threshold by integrating the joint p.d.f. of the relevant randomly distributed variables over the volume defined by the condition > * .
Regarding the BH population, the random variables entering the flux expression are: speed, mass and distance from Earth. We neglect any possible correlation between these variables. In addition, we also treat as a random variable the density of the insterstellar medium at each BH location. In both analysis we present below, we make the assumption that the gas density and BH position are not correlated. Therefore, we obtain the expected number of luminous sources sources by computing the following integral: Here (v BH ), ( ), ( ) and ( ) are the normalized p.d.fs describing respectively the BH speed, mass, distance from Earth and the interstellar medium density at its location. With { } we indicate the free parameters entering the expression for the flux, and which we consider to be fixed for all BHs. These are: the sound speed of the neutral medium , the sound speed in the ionized region , the fraction of bolometric luminosity / and, in the case of BHL accretion, the suppression factor (regarding the first of these, we remark that variations of the sound speed of the neutral medium have a negligible impact on the flux magnitude). Finally, the parameters that define the p.d.f.s for the the random variables, which we indicate with { }, should also be regarded as free parameters of the model. We discuss them in the next section.

CHARACTERIZING THE BLACK HOLE POPULATION
Mass. We assume BH masses to follow a normal distribution of mean mass = 7.8 and mass = 1.2 . This distribution was obtained from the study of X-ray binaries by Özel et al. (2010) and confirmed in an independent analysis by Farr et al. (2011). It has also been employed by previous studies on this topic such as the works of Fender et al. (2013) and Tsuna et al. (2018).
Speed. The BH velocity is given by the combination of two components: the velocity of the progenitor star, v star , and the kick the BH receives at birth due to the supernova explosion, v kick . As for the former, we are concerned with the velocities relative to the molecular gas, so we ignore the rotational component along the Galactic disk and consider exclusively the velocity dispersion . This is well measured for both nearby stars and stars in the Galactic centre.
On the other hand, little is known about the magnitude of the natal BH kicks. A study of pulsars proper motions by Hobbs et al. (2005) concluded that neutron star kicks obey a Maxwell Boltzmann distribution with = 265 km/s. Re-scaled with the average masses (Fender et al. 2013), this gives kick ≈ 50 km/s, corresponding to a mean of kick ≈ 75 km/s. We consider here as reference values kick = 50 km/s ("low kick") and kick = 100 km/s ("high kick").
The BH speed, resulting from the combination of these two independent components, is then distributed following a Maxwell Boltzmann of mean: BH = √︃ 2 star + 2 kick ,. Hence, the uncertainties in both the kick and the progenitor star's velocity dispersion are enclosed in only one parameter, BH . Distance from Earth. We will assume the BHs to be distributed uniformly in both regions under analysis and derive the distance distribution accordingly. Density of the interstellar medium. The interstellar medium is usually described as composed by three different components: the ionized component, the neutral gas, and the molecular clouds. In the PR13 picture the probability of obtaining large enough fluxes from a BH found in the first two of these components turns out to be negligible. We will therefore consider only densities associated to molecular clouds.

SEARCHING FOR BLACK HOLES IN THE LOCAL REGION
We start by applying our model to the study of the innermost 250 pc around Earth, following the work of Fender et al. (2013). In this study, the authors applied the Bondi formula to model accretion and estimated the number of visible sources through Montecarlo simulations. They considered a few combinations of the values of the mean BH speed BH and of the suppression factor , obtaining a bound on the combination of these parameters. Here we reproduce the same setup to obtain a prediction for the number of detectable Xray sources, verifying the result for the BHL model and comparing it to the predictions of the PR13 model.

Our setup
Regarding the BH population, we assume local = 3.5×10 4 objects are present in the region between 70 pc < d < 250 pc (Fender et al. 2013). We assume a uniform spatial distribution, a normal mass distribution and a MB velocity distribution of average BH , as described in the previous section. As for the interstellar medium, we assume two types of molecular clouds to be present, each type with uniform density. These are: a) warm clouds, of number density 10 2 cm −3 , sound speed 0.9 km/s, filling factor 5 × 10 −2 , and b) cold clouds, of number density 10 3 cm −3 , sound speed 0.6 km/s and filling factor 5 × 10 −3 (Fender et al. 2013). This combination results in an average density avg = 10 −3 . The integral in Eq. 15 becomes: Here clouds indicates the number of BHs we expect to find in each type of cloud: where clouds is the associated filling factor. We obtain clouds ≈ 175 in the cold clouds and clouds ≈ 1750 in the warm ones. Finally, the bolometric fraction of luminosity in the X-ray band is set to 0.3.

Results
Let us first consider the X-ray flux associated to a generic isolated BH located in the local region and accreting gas from a molecular cloud as a function of the BH speed. In figure 2 we show this flux for both the PR13 accretion scenario and the conventional modeling based on the BHL formalism, the latter suppressed by a factor = 0.01. The two colored bands correspond to the two densities considered for molecular clouds. The width of the band is obtained considering distances between 70 and 250 pc and masses between 5 and 11 solar masses. We show for reference the detection threshold associated to the 4th INTEGRAL IBIS/ISGRI soft gamma-ray survey catalogue (17-100 keV; Bird et al. 2009) , which is approximately 10 −11 erg/cm 2 /s. We want to emphasize once again the distinct behaviour associated to the two accretion scenarios. In the Bondi picture, the slower sources are very luminous. This is particularly relevant in this context, since velocities in the local region are expected to be on the low end of the range shown in this figure: the unsuppressed BHL scenario predicts a huge number of sources (see figure 3) and can easily overshoot the bound. As a consequence of the introduction of a suppression factor, only very slow sources are expected to be visible. This is no longer true when the PR13 scenario is considered, since it naturally predicts a flux suppression for the slower sources. Instead, the population of BHs emerging above threshold and showing up in the X-ray sky will present relatively high speeds (around 50 km/s). Furthermore, we can notice how, based on the Bondi picture, we expect to detect BHs in both types of clouds. According to the PR13 model, however, only BHs located in the denser clouds can be detected. Following the procedure described in the previous section, we can now integrate the distributions that characterize the BH population to obtain the number of sources visible in the X-ray sky with the IBIS/IGRI survey. In figure 3 we show the dependence of number of sources on the average speed BH . This prediction is obtained for three accretion scenarios: a) the PR13 model; b) the BHL model with suppression factor = 0.01 c) the BHL model with perfect efficiency = 1. We indicate with orange and green arrows the reference values given by the "low kick" the "high kick" scenarios, considering a velocity dispersion of the progenitor stars of 15 km/s (corresponding to an average speed of ≈ 25 km/s). For comparison, the average speeds corresponding to the "no kick" and "high kick" ( kick = 75 km/s) scenarios considered in Fender et al. (2013) are also shown and indicated with A and B, respectively.
Our results can be summarized as follows: • As far as the BHL scenario is concerned, our findings are in agreement with the ones obtained by the results reported in Fender et al. (2013): i.e., a few tens of visible sources in scenario A (no kick, = 0.01 ). Larger values of are excluded unless very high kicks are present. This is, once again, due to the rapid increase of the BHL rate at low speeds.
• In the PR13 scenario, the introduction of a suppression factor is not necessary. While being compatible with the experimental bound, our model nevertheless predicts a significant number of bright sources. This prediction corresponds to a population of ∼ 40-100 accreting isolated black holes in the existing catalogues, taking as reference the "no kick" and "high kick" scenarios . We will discuss in more detail the consequences of this result in the next Sections, and compare this prediction with the one associated to the Central Molecular Zone.

SEARCHING FOR BLACK HOLES IN THE CENTRAL MOLECULAR ZONE
In this section we turn our attention to the inner part of the Galactic bulge. This region is a particularly promising target for BH searches because of the presence of a large reservoir of molecular gas called central molecular zone (CMZ), a cloud complex with asymmetric shape that extends up to 200 pc away from the Galactic centre. The CMZ has been the object of extensive multi-wavelength observational campaigns over the years, aimed at characterizing its structure and physical properties, including star formation rate (see for instance Kruĳssen et al. (2019) and references therein for a recent analysis). The abundance of gas in the form of giant molecular clouds, and the location close to the gravity centre of the Galaxy clearly makes it the ideal target for our study.

Our setup
To model the CMZ, we employ a simplified version of the model by Ferrière et al. (2007). We describe it as a cylindrical region of half height 15 pc and radius 160 pc. The number of BHs contained in this region is of course very uncertain, but we can make a naive order-of-magnitude estimate, as follows. We assume the BHs to be generated following the distribution of stars in the Galaxy. Here we are interested in the bulge component, which accounts for around 15% of the total and can be modelled with a spherical exponential with scale radius bulge = 120 pc (Sofue 2013). Integrating the bulge spherical exponential distribution over the CMZ volume gives us the fraction of ABHs born in this region: around 2% of the bulge component. Assuming a total of 10 8 BHs present in the whole Galaxy, this corresponds to 3 × 10 5 BHs. We must however take into account that, due to the large initial natal kicks, the initial spatial distribution is modified (Tsuna et al. 2018). We performed a simulation of the evolution in the Galactic potential (Irrgang et al. (2013), model II) of 1000 BHs. These are initially distributed uniformly in the CMZ region with an average speed of 130 km/s and given an average natal kick of 75 km/s. The simulation shows that only ≈ 25% of the BHs remain in the region: we observe in particular a spreading in the direction perpendicular to the Galactic plane. On the other hand, we do expect some BHs born outside of the CMZ to enter the region under the effect of the gravitational potential. However, these objects cross the CMZ close to the periastron of their orbit, with very large velocities and hence suppressed accretion rates. We can therefore neglect the latter effect. However the former is significant and we take it into account. We thus update our naive estimate to CMZ ≈ 7.5 × 10 4 . The molecular clouds in this region are known to be denser than the Galactic average. We assume an average number density of molecular hydrogen of avg = 150 cm −3 (Ferrière et al. 2007) and consider two types of clouds: warm, less dense clouds, and cold, denser ones. The warm clouds are of secondary importance given the large distance from Earth (see figure 4). We set their density to = 10 2.5 cm −3 and consider a filling factor of 14% (Ferrière et al. 2007). Regarding the cold clouds, we choose to adopt a more refined modelling based on a power law density distribution between min and max : Based on Ferrière et al. (2007), we set min = 10 3.5 cm −3 and max = 10 5 cm −3 , while we treat as a free parameter. The filling factor for the cold clouds is then obtained by requiring avg = 150 cm −3 . For a reference value of = 2.4 (He et al. 2019(He et al. , 2020 we obtain a filling factor of ≈ 1% and an average cloud density of ≈ 8 × 10 3 .
The temperature of the clouds enters the accretion rate through the sound speed . As we mentioned before, this variable has little impact on the predictions and we set it to 1 km/s.
Regarding the speed distribution, the average speed BH is We show the X-ray flux associated to a BH orbiting in the Central Molecular Zone molecular cloud complex, for different values of the relative velocity with respect to the gas cloud, and the gas density in the cloud. We consider the Park-Ricotti and the Bondi-Hoyle-Littleton formalism. As a reference, we show the detection threshold associated to the NuSTAR Galactic Center survey (Hong et al. 2016). treated as a free parameter. Nevertheless, we can make some estimates of reasonable values. The most recent observations suggest a velocity dispersion of the central bulge stars of average bulge ≈ 130 km/s (Sanders et al. 2019, Brown et al. 2018. This gives, for reference, = 140 km/s in the "low kick" scenario and = 160 km/s in the "high kick" scenario.
Finally, we approximate the distance distribution with a delta peaked at 8.3 kpc and we assume the normal mass distribution discussed in Section 5.
The integral to be computed becomes: The expected number of black holes in the clouds is obtained as: where clouds is the associated filling factor. In summary, the set of relevant parameters that play a key role in our calculation are: 1) the average BH speed, BH 2) the slope of the power law density distribution, 3) the sound speed within the ionized region around the BH, in ; 4) the fraction of bolometric luminosity irradiated in the hard X-ray band x / ; 5) the number of black holes in the CMZ, CMZ .

Results
We start by considering the X-ray flux associated to a generic isolated BH located in the Central Molecular Zone and accreting gas from a molecular cloud as a function of the speed and of the molecular cloud density. We show this observable in figure 4 for the two accretion scenarios discussed above. As a reference, we show the detection threshold associated to the NuSTAR Galactic Center survey in the 3-40 keV band Hong et al. (2016). The same trends highlighted in Sec. 5.2 can be noticed in this plot: In particular, we point out that, if the cloud density is high enough, there is a wide range of BH speed associated to emission above threshold within the PR13 formalism. On the other hand, no sources are expected to be detected in the lower density clouds.
With the X-ray flux for a CMZ source at hand, we apply the procedure described in the previous Sections, and compute the number of X-ray sources associated to accreting BHs in this region. The results are shown in figure 5. In this panel we show the impact of the different free parameters on this key observable. In particular, we show in panel (a) the number of X-ray sources in the CMZ region as a function of the average speed of the BH population adopting both the BHL and the PR13 accretion models. We show the result for different choices of the sound speed in the ionized region, and for different choices of the parameter in the BHL case. We can notice how, in this setting, the PR13 scenario predicts as many sources as the = 1 BHL, due to the fact that high speeds are prevalent among the BH population of the CMZ. However, we have seen that the the un-suppressed BHL scenario is excluded by previous studies on nearby compact objects and complementary studies focused on AGN populations as mentioned in Sec. 2.1. The comparison should therefore be made between the predictions of the PR13 model and the suppressed BHL model, also shown in figure for = 0.01. Then, the PR13 model predicts a significantly greater number than those obtained assuming BHL accretion. (See also the work of Tsuna et al. (2018), in which O(1) sources were predicted in the Galactic centre using BHL accretion). In panel (b) we focus on the dependence of the same observable with respect to the fraction of the bolometric luminosity that is radiated in the X-ray band of interest; in panel (c) we consider the differential contribution of the clumps of different density, again for different reference values of the ionized sound speed in the BH vicinity . We can conclude from the results visualized in figure 5 that, despite the relevant uncertainties associated to the modelling of such a complex setup (and despite the existence of a threshold effect due to the peak in the accretion rate), the prediction of few tens of sources is solid with respect to these uncertainties. In particular, we notice from panel (a) that a number of bright X-ray sources comprised between 10 and 20 is expected for our reference value of and for = 150 km/s. As far as the speed distribution is concerned, we remark that the two reference values we have chosen to bracket the uncertainty both lead to similar predictions, while large deviations from this range would require to assume unrealistically high speeds. The same line of thought applies to the slope associated to the clump density distribution. From panel (b) we also notice a lower number of sources may only arise if we were to assume either a very high temperature in the ionized region, or else a very low fraction of bolometric luminosity radiated in the X-ray band. We remark,however, that the prediction scales linearly with the expected number of BHs in the region, CMZ , which carries a large uncertainty.
In summary, the PR13 model (together with our rather conservative assumptions regarding the X-ray emission) leads to the prediction of a significant number of X-ray sources in the CMZ region associated to isolated BHs accreting from molecular clouds. The total number of sources is of the same order of magnitude as the one predicted in the local region and discussed in the previous Section. However, while the previous Section was dealing with a full-sky analysis, in this case we are focusing on a specific region of interest (ROI) with small angular extent, which is ideal for multiwavelength observational campaigns. We recall that NuSTAR has detected 70 hard X-ray sources sources in this region interest. The nature of most of the sources in the catalogue in not precisely known,   although a significant population of cataclysmic variables and Xray binaries is expected. Our finding seems to suggest the need of a careful data analysis, in order to identify the possible presence of a relevant population of isolated accreting black holes in the already existing data. Following this line of thought, and in the prospect of a discovery, a multi-wavelength analysis is compelling, so we will turn our attention to the radio domain in the next Section.

MULTI-WAVELENGTH PROSPECTS FOR THE SQUARE KILOMETER ARRAY
In the previous section we have shown that, in a wide portion of the parameter space associated to our problem, a large number of bright X-ray sources are expected in the Galactic Center region. However, in order to pinpoint a source as an accreting black hole, a careful multi-wavelength study has to be performed. The GHz radio band is particularly interesting in this context. As in previous works, (Gaggero et  scaling between X-ray and radio flux (Plotkin et al. 2012b) discussed in detail in Section 3 (eq. 14) is employed to obtain the radio flux. A remarkable increase in the sensitivity is expected in the radio domain over the coming decade, thanks to the development of the Square Kilometer Array (SKA) project. This experiment has a huge potential towards shedding light on key problems of fundamental physics, cosmology and astrophysics (Weltman et al. 2020). Here we focus on the discovery potential of the population of astrophysical black holes in the CMZ at 1.4 GHz band, and provide estimates of the number of potentially detectable sources by the SKA1-MID facility. Assuming gain = 15 K/Jy, receiver temperature rx = 25 K, sky temperature towards the Galactic Center sky = 70 K, and bandwidth = 770 MHz (as in Calore et al. 2016), we obtain an instrumental detection sensibility of 2.7 Jy for one-hour exposure. In the following, we assume an optimistic 1000 h exposure time and consistently adopt a 85 nJy as potential detection threshold.
In figure 7 we show the radio flux associated to a generic BH in the CMZ cloud, and compare it with the prospective SKA sensitivity and with the threshold of the VLA catalogue (Lazio & Cordes 2008). Focusing on the PR13 scenario, we can conclude from this plot that, while we do would not expect to detect any isolated BH give the VLA sensitivity , SKA would be able to unveil a huge population of isolated BH in the Galactic centre. In fact, most of the BHs accreting from the cold clouds would emit above the SKA threshold.
We show our prediction for the number of radio sources detectable with SKA in figure 8 and figure 9, obtained with the same setup described in section 7.1. In particular, in figure 8 we show the number of radio sources detectable by SKA as a function of the BH speed, for different choices of the parameters discussed in the previous Section. We can notice how our model leads to predict thousands of visible sources the two reference values we have chosen to bracket the uncertainty on the speed distribution.
Let us now widen our perspective and promote CMZ , the total number of BHs in the CMZ to a free parameter. In figure 9 we show the number of bright radio sources in the ( , CMZ ) space. In this parameter space, we may use the comparison with NuSTAR observations discussed in the previous Section to obtain a bound: in particular, we can identify an upper limit on the number of black holes present in the CMZ, by requiring not to overshoot the number of sources observed by that experiment. We find that, in the allowed portion of parameter space, a very large number of radio point sources will be detectable by SKA. In summary, the radio searches for isolated BHs in the CMZ constitute a promising science case for this experiment.

DISCUSSION
Comments on the parameter space. We have concluded that our predictions for the CMZ are stable in a very wide region of the parameter space of the model. However, we should note that we made a number of assumptions regarding the distribution of black holes in the CMZ. On the one hand, as we already discussed, our estimate of the total number of BH in the CMZ is obtained with a simple model for the initial distribution and subsequent orbit evolution. Furthermore, we have assumed that both the BHs and the clouds are uniformly distributed in the CMZ. The presence of a positive or negative correlation between the two distributions could in principle impact the results. We further assumed a specific model for the cloud density distribution. These simplifications eventually affect the number of BHs in the clouds, and therefore the final predictions on the number of sources, which scale linearly with this quantity.
The accretion scenario. We have shown in Section 6 that, based on local observations, the PR13 accretion model does not require the introduction of a suppression factor . However, since the PR13 model relies on the BHL model to describe accretion within the ionized region (see Sec. 2), one may wonder whether the BHL suppression factor should be introduced at that stage. With respect to this, we should take into account that this suppression factor is usually associated to a very strong outflow of matter from within the Bondi radius, associated to a fraction of 1 − of the total matter initially accreted. Such a large outflow is not observed in the simulations of Park & Ricotti (2013), where most matter crossing the Bondi radius within the ionized region ultimately reaches the BH. We conclude from this that introducing a suppression factor in the expression would not be physically justified. We remark that the simulations considered here represent the state of the art in the field. However, some physical ingredients are still not captured: we mention in particular the role of magnetic fields. Future works will assess the potential impact of a full magneto-hydrodynamical treatment in this context.
The emission mechanism and spectrum. In Sec. 7.1 we also assume that ∼ 30% of the total bolometric luminosity of all our isolated BHs falls within the observable X-ray band-X = 0.3 bol . This assumption was based on previous estimates, for example those presented by Fender et al. (2013). However, it is worth assessing the accuracy and precision of this assumption. We assumed a flat radio spectrum ( ∝ 0 ), and then an additional inverted power law in the X-ray band. The assumption that 30% of the observed bolometric luminosity falls in the X-ray band depends on the spectral index of this additional power law, , the high-energy cut off of the power law, cut , the absorbing hydrogen column density along the line of sight, H , and the observing energy band (which we conservatively adopted as 3-40 keV for the NuSTAR X-ray telescope). Inefficiently radiating Galactic BH-XRBs are dominated by power law emission with a typical range for the spectral index given by 0.6 ≤ ≤ 1.0 (or a photon index of 1.6 ≤ Γ ≤ 2 in common X-ray astronomy parlance). However, the power law spectra of Galactic BH-XRBs can become softer in intermediate spectral states, or very quiescent states, so we can extend this out to an upper bound of ≤ 1.4 or Γ ≤ 2.4. Adopting a rough estimate of the jet break frequency of ∼ 10 14 Hz (infrared frequencies; Russell et al. 2013), a typical highenergy cut off in the X-ray of 100 keV, and a high hydrogen column density of 10 23 cm −2 , we calculate that X ∼ (0.24-0.35) bol . The corresponding unabsorbed range of fractional luminosities is then 0.27-0.46 bol . These predictions, however, depend strongly on the cutoff energy. For example, just reducing the cutoff energy from 100 keV to 80 keV can increase the fractional luminosity to ∼ 60%-70%. Therefore, our assumption that X ∼ 0.3 bol is rather conservative, even if we assume strong absorption along the line of sight. It is also based upon a relatively hard power law spectrum in the X-ray band. Nevertheless, we have shown in Section 7, that our predictions are solid with respect to variations of this parameter, even for values as large as 70%.
For the radio flux calculations through the Fundamental Plane of Black Hole Activity (FP; Merloni et al. (2003); Falcke et al. (2004); Plotkin et al. (2012b)) we had to assume the presence of a jet emitting in the GHz radio band, following Fender et al. (2013); Gaggero et al. (2017); Manshanden et al. (2019). This crucial assumption is motivated by results of Fender (2001), which established that all highly sub-Eddington accreting BHs in X-ray binary systems are accompanied by such radio emission.
Further emission could come from other types of outflow that are less collimated or relativistic. Such alternative emission could constitute a lower limit on the expected radio emission, which would be of particular interest in the absence of a jet. Tsuna & Kawanaka (2019) assumed a spherically symmetric outflow based on losing a fraction (1-) of the accreting matter before reaching the BH, where is the suppression factor. Assuming a power-law radial dependence of the accretion rate, they used the escape velocity to conservatively estimate the power associated to such an outflow. It should be noted that this radial power-law assumption (with slope larger than 1) implies that most matter would be lost close to the BH, which requires a powerful outflow. The shock from the collision of this outflow with the ISM could accelerate non-thermal electrons through diffusive shock acceleration and subsequently produce radio-synchrotron emission. The precise emission profile will be highly dependent on the fraction of energy going into the electrons and the magnetic field.
However, as highlighted above, such a relevant outflow is not observed in the simulations presented in Park & Ricotti (2013) and hence this scenario is hard to unify with the PR13 accretion model. Instead, it would be interesting to further investigate what type of emission could be associated to the shock on the boundary of the cometary region itself, which would likely require additional detailed simulations.
Connection with PBH searches. We now want to widen the perspective and mention that the search for isolated astrophysical black holes (ABHs) can be considered a relevant problem per se, and also crucial as a background in the context of the quest for primordial black holes (PBHs) in the Galactic environment. The potential of this observational channel in constraining the abundance of PBHs in the inner part of the Galaxy, with particular focus on the CMZ region, was first pointed out in Gaggero et al. (2017), within the framework of the Bondi-Hoyle-Littleton formalism. In Manshanden et al. (2019) the impact of the Park-Ricotti formalism was assessed in detail. The main result of that work is a strong upper bound on the abundance of PBHs (in the context of a non-clustered distribution of PBHs with monocromatic, log-normal and powerlaw mass functions). The upper bound was set at the level of 10 −3 of the Dark Matter in the form of PBHs, which would correspond to roughly ∼ 10 8 objects in the Milky Way for a 10 M reference mass. Hence, the astrophysical population (which is expected to include the same order of magnitude of compact objects) is expected to be an irriducible background that could play the role of an astrophysical floor associated to the quest for a subdominant population of PBHs of tens of solar masses.
Outlook. We have shown that a scenario based on the PR13 accretion model and our best estimate of the number of BHs in the CMZ naturally implies a large number of detectable X-ray and radio sources, for any reasonable choice of the other free parameters involved in the problem.
This relevant result is highly suggestive that a clear identification of a population of isolated, accreting BHs in the Galactic Center region -by means of the current or forthcoming generation of X-ray experiments -may be around the corner.

SUMMARY
In this paper we presented a comprehensive study of the constraints and prospects of detection of a population of isolated astrophysical black holes in the solar vicinity and near the Galactic Centre. We have adopted the Park-Ricotti (PR13) accretion model that takes into account radiation feedback and is backed up by hydrodynamic numerical simulations, and compared the results obtained within this scenario with the ones obtained by exploiting the well-known Bondi-Hoyle-Littleton formalism. We found that, within the PR13 model, a large number (∼ 40-100) of bright X-ray sources associated to isolated accreting black holes is expected in the solar vicinity ( < 250 pc), and a significant fraction of the sources listed in current catalogues could be associated to such objects. We performed an extended parametric study about the number of detectable sources in the vicinity of the Galactic Centre, with particular focus to the promising Central Molecular Zone region. In this region of interest we found that, despite the large uncertainty associated to the free parameters involved in the calculation, the prediction of O (10) Xray sources above the detection threshold associated to the NuSTAR catalogue is solid for any reasonable choice of the parameters under scrutiny. In particular, a number of bright X-ray sources comprised between 10 and 20 is expected for our reference value of the ionized sound speed = 30 km/s and an average black hole speed BH = 150 km/s. A multi-wavelength analysis will be needed to clearly identify such a population. In this context, we pointed out promising prospects of detection for future generation experiments in the radio band, with particular reference to the Square Kilometre Array.