Searching for Primordial Black Holes in the radio and X-ray sky
We model the accretion of gas onto a population of massive primordial black holes in the Milky Way, and compare the predicted radio and X-ray emission with observational data. We show that under conservative assumptions on the accretion process, the possibility that primordial black holes can account for all of the dark matter in the Milky Way is excluded at by a comparison with a VLA radio catalog at GHz, and at by a comparison with a Chandra X-ray catalog ( keV). We argue that this method can be used to identify such a population of primordial black holes with more sensitive future radio and X-ray surveys.
Introduction: The first direct detection of a gravitational wave signal, announced by the LIGO collaboration earlier this year Abbott:2016blz demonstrated the existence of black holes (BHs), prompting the suggestion Bird:2016dcv; Clesse:2016vqa that these objects are primordial black holes (PBHs) that may account for all of the dark matter (DM) Jungman:1995df; Bertone:2004pz; Bertone:2010at in the Universe. The connection between PBHs and DM has been extensively studied in the past (see e.g. Ivanov:1994pa; Khlopov:2008qy; Carr:2009jm; Blais:2002nd; Afshordi:2003zb; Frampton:2010sw), and a number of constraints exist on the cosmic abundance of PBHs over a very wide mass range (see the discussion below, and e.g. Ref. Carr:2016drx for a recent review).
In this Letter, we consider for the first time, in the context of PBH searches, the X-ray and radio emission from the Galactic Ridge region produced by the accretion of interstellar gas onto a population of PBHs in the Milky Way. Given current estimates of the bulge mass Portail2015, if PBHs constitute all of the DM, there should be such objects within kpc from the Galactic center (GC). Since the inner part of the bulge contains high gas densities Ferriere2007, a significant fraction would inevitably form an accretion disk and emit a broad-band spectrum of radiation. We show (fig. 1) that radio and X-ray data in the Galactic Ridge region rule out, at 5 and 40 respectively, the possibility that PBHs constitute all of the DM in the Galaxy, even under conservative assumptions on the physics of accretion.
Our limits arise from a realistic modeling of the accretion process, based on the observational evidence for inefficient accretion in the Milky Way today Perna:2003; Pellegrini:2005, and corroborate, with a completely independent approach, the exclusion of massive PBHs as DM candidates.
Accretion on black holes: We should expect the accretion rates, , of a Galactic population of PBHs accreting from interstellar gas to be well below the Eddington limit . Even under the unrealistic assumption of Bondi-Hoyle-Lyttleton accretion Hoyle1939; Bondi1944, and typical velocities as low as km/s, the accretion rate would definitely be sub-Eddington: .
BHs accreting at are radiatively inefficient, such that the luminosity scales non-linearly with heinz2003. The prevailing physical pictures adopted to explain the weak emission properties are advection-dominated accretion in which the gas cooling timescales greatly exceed the dynamical timescales Narayan:1994, and mass loss from the disk or internal convective flows, such that the accretion rate itself has decreased once gas reaches the inner edge of the disk blandford:1999; quataert2000. It is likely that both mechanisms are at play, a view supported by both radio and X-ray constraints on the gas density around Sgr A*, the supermassive BH at the center of the Galaxy, the least luminous accreting BH observed to date (in Eddingtion units), and thus a well-studied source from the point of view of weak accretion physics Bower2003; Marrone2007; wang2013. We compute the accretion rates and the radiative efficiencies of a Galactic population of PBHs in the low-efficiency limit, following the formalism presented in Maccarone2005MNRAS; Fender:2013. We take into account the findings of previous studies regarding accretion of interstellar gas onto isolated black holes Fujita1998; Armitage1999; Agol2002.
We model the radiative efficiency , defined by the relation for the bolometric luminosity , as for (if we were to assume instead efficient accretion above the critical rate, , then we would have a constant ). As already discussed, all our sources fall below this critical accretion rate, such that they are all inefficient accretors: this means the luminosity scales non-linearly with accretion rate, .
We parameterize the accretion rate as , such that
where is the gravitational constant, is the velocity of the BH, and is the sound speed of the accreted gas, which is below 1 km/s in cold, dense environments.
An important element that needs consideration is the temperature of the accreted gas due to radiative pre-heating Maccarone2005MNRAS. Photoionising radiation will lead to an ionisation bubble surrounding the source, known as the Strömgren sphere Stromgren1939, with a characteristic radius, . In the following, we assume that the gas around the BH is fully ionized – and therefore, we set km/s – if the timescale for the ionization of the Strömgren sphere is shorter than the timescale associated with the incoming flux of fresh, unprocessed material.
Regarding , we choose a reference value of . Given the degeneracy between and the angular momentum and temperature of the accreted gas, this value is consistent with isolated neutron star population estimates and studies of active Galactic nuclei accretion Perna:2003; Pellegrini:2005; wang2013
This prescription is the same as that adopted by Fender:2013; however, we consider , and rescale the value of used in that work across the full 10–100 mass range.
We convert bolometric luminosity to X-ray luminosity via the approximate factor following Fender:2013.
Motivated by the results presented in Fender:2001 and by the discussion in Maccarone2005MNRAS; Fender:2013,
we assume the presence of a jet – thus requiring a system with a surplus of angular momentum, or a dynamically important magnetic field combined with a spinning black hole –
emitting radio waves in the GHz domain with an optically thick, almost flat spectrum, whilst the
X-ray emission is non-thermally dominated, originating from optically thin regions closer to the
BH. In order to convert the X-ray luminosity into a GHz radio flux, we adopt the universal empirical
relation discussed e.g. in Plotkin:2012, also known as the fundamental plane (FP), which
applies for a remarkably large class of compact objects of different masses, from X-ray binary systems to active Galactic nuclei. We calculate
the X-ray luminosity in the 2–10 keV band in accordance with the FP, assuming a hard power-law X-ray spectrum with photon index ,
and a typical range for hard state X-ray binaries of 1.6–2.0 (see Hong:2016).
We extrapolate this power-law spectrum into the 0.5–8 keV and 10–40 keV bands in order to also make comparisons with Chandra and NuSTAR catalogs.
We then use the FP relation to calculate the 5 GHz
radio flux from the 2–10 keV X-ray flux and assume a flat radio spectrum, such that , allowing direct comparison with the 1.4 GHz source catalog from a VLA survey of the GC region.
Primordial black hole population: In order to derive a bound from X-ray and radio data, we set up a Monte Carlo simulation for each PBH mass, assuming a delta mass function.
We populate the Galaxy with PBHs following the Navarro-Frenk-White (NFW) distribution Navarro:1995iw (other more conservative choices are discussed below). We implement the accurate 3D distribution of molecular, atomic, ionized gas in the inner bulge presented in (Ferriere2007); that distribution includes a detailed model of the 3D structure of the Central Molecular Zone (CMZ), a 300 pc wide region characterized by large molecular gas density and centered on the GC, i.e. in the region where the largest density of PBHs is expected.
For each PBH, the velocity is drawn randomly from a Maxwell-Boltzmann distribution. The characteristic velocity of the distribution is position-dependent. The velocity distribution at a given radius is a crucial ingredient, because the accretion rate scales as , eq. (1). In order to derive such a distribution, we consider the recent state-of-the-art model for the mass distribution in the Milky Way described in McMillan2016, where 6 axis-symmetric components are taken into account (bulge, DM halo, thin and thick stellar discs, and HI and molecular gas discs). We then assume that the velocity distribution at a distance from the GC is a Maxwell-Boltzmann with . Under the assumption of isotropic orbits111We verified that, in the high-resolution Aquarius N-body simulations, the anisotropy parameter is consistent with 0 in the whole range of radii we are interested in, therefore the assumption of isotropic orbits is solid., an exact computation of the phase-space density could be performed by means of the Eddington formalism Eddington:1915, as done e.g. in Fornasa:2014. We checked that our simple approach is equivalent in the low-velocity tail, up to km/s. 222M. Fornasa, private communication. Since our results depend only on PBHs with velocities km/s (see below), we can safely neglect the high-velocity tail and adopt the simple formalism described above.
Given the mass, position and velocity of each PBH (and the gas density), we compute accretion rate, X-ray, and radio emission adopting the prescriptions discussed in the previous section.
Radio BH candidates: The GHz source catalog from a VLA survey of the GC region Lazio2008 contains sources in a region centered on the GC. The minimum detectable flux for this catalog is mJy.
In order to compare our predictions to the observations, we carry out a data analysis on the VLA catalog and check if there can be any BH candidate among the detected sources.
If any of these sources are accreting BHs, their X-ray and radio emissions should be co-located. We therefore compare the radio catalog with the X-ray point source catalog from Muno2009, which contains sources detected by Chandra in the keV band in a band centered on the GC, and search for all sources in both catalogs that have positions within of each other.333This is a very conservative separation. The positional accuracy of Chandra is . For the VLA, the positional accuracy is typically a small fraction of the synthesized beam, for the survey in Lazio2008 , taken in A configuration. A separation of is chosen in Lazio2008 to search for positional coincidences in other radio catalogs; we therefore also adopt as the maximum allowed separation.
We find sources in both the X-ray and radio catalogs within of each other. If we assume that these sources are accreting BHs, then their X-ray and radio fluxes should lie on the FP, as explained above. So, we use the FP ( considering masses from 10 to 100 ) to predict the X-ray flux from the radio flux of each of these objects ( in the very conservative case, if we exclude likely foreground sources).
We find that the predicted X-ray fluxes are substantially larger ( orders of magnitude) than the flux reported in the catalog from Muno2009 in the whole mass range we consider. We therefore conclude that none of the 24 (or 9 likely Galactic) VLA sources with overlapping positions lie on the FP, and therefore, given the assumptions described above regarding the presence of a jet, we have no BH candidate in our sample.
X-ray BH candidates: Hard X-ray emission ( keV) suffers from far less Galactic absorption than soft X-ray emission and is therefore a good band to search for emission from accreting BHs.
We consider sources in the Chandra catalog Muno2009 in the keV band, and those detected by NuSTAR in the keV band Harrison:2013. For Chandra (NuSTAR), we consider a small region-of-interest (ROI) including the high-density region of the Galactic Ridge: (). There are likely Galactic X-ray sources in the Chandra catalog above a flux threshold of ph cm s444”Likely Galactic” sources are defined in Muno2009 based on their hardness ratios. The exposure across the Chandra survey region is variable and the flux threshold used here is a compromise between maximizing the ROI and the completeness, per Muno2009., and NuSTAR sources. Since in all cases the corresponding radio flux predicted with the FP would be orders of magnitude below the detection threshold of the VLA survey in Lazio2008, we cannot draw any conclusions on the nature of these X-ray sources. Therefore, we consider all of them in our analysis as potential BH candidates (we only remove of the detected NuSTAR sources that are thought to be cataclysmic variables Hong:2016).
Results: The main result of the Letter is presented in fig. 1. We display the 2, 3, and 5 constraints on the DM fraction as a function of the PBH mass.
The upper limits are derived as follows. We perform Monte Carlo simulations for reference values of the mass in the interval, assuming a DM fraction . We determine the mean and standard deviation of the distributions of the predicted number of PBHs with radio fluxes above the VLA threshold and with X-ray fluxes exceeding the Chandra (NuSTAR) threshold, in the corresponding ROIs. We verify that the number of bright PBHs is compatible with Poisson statistic and the average predicted number scales linearly with . We derive the radio and X-ray bounds by comparing the number of predicted PBHs with the number of BH candidates derived from the analysis of radio and X-ray catalogs described in the previous section. For the X-ray bound, we show the result obtained with the more sensitive Chandra catalog. The NuSTAR bound is slightly weaker: It allows us to exclude at values of as low as (for ).
In fig. 2, we show the PBHs detectable by VLA at GHz assuming a PBH mass of 30 and DM fraction equal to 1, for one specific Monte Carlo realization. This scenario predicts, on average, 40 6 sources above the VLA flux threshold for and, thus, it is excluded by more than 5 from radio observations. However, it is important to understand where the constraining power comes from: The PBHs above the detection threshold, and thus the ones with the larger X-ray flux, lie in the very inner region of the Galaxy where the column gas density is the highest and show very small velocities, in the range 0.3 3 km/s. Therefore, the constraints arise from the very low velocity tail of the distribution and from regions correlated with very high column densities, e.g. CMZ, as already mentioned above.
Discussion and conclusions: In this Letter we derive new, strong constraints on the hypothesis that PBHs comprise all of the DM in the Universe. In particular, we find that PBHs with , that could be responsible for the gravitational waves detected by LIGO, contribute less than 20% to the whole DM density.
In the mass window , our constraints are competitive with (and even stronger than) those arising from the study of microlensing events with the MACHO project Alcock:2001 (for ) and with those from halo wide binaries 2009MNRAS.396L..11Q; 2014ApJ...790..159M (for ). For , they are also comparable or stronger than the constraints from the survival of central star clusters in faint dwarf galaxies, in particular in Eridanus II Brandt:2016aco; Li:2016utv. Even more stringent constraints arise in principle from the analysis of the Cosmic Microwave Background (CMB) Ricotti:2007au. However, those arising from the analysis of spectral distortions (based on FIRAS data) turned out to be much weaker than originally thought ClesseBellido2016, while the ones based on the study of CMB anisotropies (see also the recent results by Chen2016), are based on assumptions on the accretion of gas on PBHs in the early Universe that are still under debate, as the modeling of the accretion process is based on theoretical arguments, and not directly supported by observations as in our case (see also the discussion in ClesseBellido2016).
In contrast with Ricotti:2007au, in fact, we adopt a very conservative prescription, compatible with current astronomical observations, for both the accretion rate and the radiative efficiency, setting the ratio of the actual accretion rate to the Bondi rate, , equal to 0.02. We remark that probably follows some distribution and is also likely degenerate with and —future studies are required to disentangle these. Moreover, we exploit for the first time in this context the empirical FP relation between radio and X-ray emission, which has been observed on a wide class of sources in a large mass range, from X-ray binaries to active Galactic nuclei. By adopting such a relation, we are able to predict the expected radio and X-ray luminosities of a population of PBHs in the Galaxy compatible with the DM phase-space distribution, as well as to look for BHs candidates in radio and X-ray catalogs. We set upper limits on the DM PBH fraction using both radio (VLA) and X-ray (Chandra, NuSTAR) point-like source catalogs, by comparing the number of expected PBHs above observational thresholds and the observed number of BH candidates in a very narrow region about the GC.
These bounds are robust with respect to the modeling of the full velocity distribution, since the predicted number of bright PBHs only depends on the very low-velocity tail ( 10 km/s) where we checked the agreement among different numerical/analytical methods. Moreover, our limits are independent of the details of the gas distribution (we checked that the bound is still present even with a naive modeling of the CMZ as a sphere with uniform density compatible with the mass constraints provided in Ferriere2007). They are also not affected significantly by a shallower DM profile as proposed e.g. in Calore2015; however, assuming an even flatter profile like the Burkert one (an extremely conservative assumption for our Galaxy), the bound is present only under the assumption of Bondi accretion.
We recall that our limits hold for a narrow mass function; a detailed study of the impact of different mass distributions is beyond the scope of the present paper and postponed to a future work.
Although our radio and X-ray bounds vanish for , future instruments will be able to verify better the accretion model as well as the PBH DM fraction. In particular, given the significant increase in sensitivity of future radio telescopes, we expect an important part of the yet-allowed parameter space to be probed by upcoming facilities such as MeerKAT and, later, SKA. Using the radiometer equation (1984bens.work..234D), the minimum (1) detectable radio flux is . For SKA1-MID (1.4 GHz), we assume gain G = 15 K/Jy, receiver temperature = 25 K, sky temperature towards the GC = 70 K, and bandwidth = 770 MHz Calore:2015bsx.
For one-hour exposure, the instrumental detection sensitivity of SKA1-MID turns out to be Jy (significantly above the source confusion limit), which would give O() detectable PBHs for our reference value (= 1, and ).
SKA will therefore be able to either place a very strong constraint in absence of BH candidate detection, or detect a subdominant population of PBHs (although the expected population of astrophysical BHs becomes comparable with the primordial one for DM fractions lower than ). With an even longer exposure ( h, 85 nJy sensitivity), achievable for dedicated deep field continuum observations, such strong constraints can be placed by SKA even under the assumption of extremely low values of , .
Interestingly, our procedure can also be applied in order to extend work on searches for astrophysical BHs in the Galaxy Agol2002; Maccarone2005MNRAS; Fender:2013, adopting the realistic spatial and velocity distributions expected for those objects. Our formalism has the potential to characterize this guaranteed population of objects in future analyses.
Acknowledgments: We are indebted to M. Fornasa for verifying of the phase-space velocity at low radii with the Eddington formalism. We thank P. D. Serpico for his useful comments on the manuscript, and D. Baumann, P. Crumley, B. Freivogel, A. King, A. Urbano, J. Vink, and C. Weniger for useful discussions. Finally, we thank the anonymous reviewers for their insightful comments and suggestions.