{'fraction_non_alphanumeric': 0.034783714608206284, 'fraction_numerical': 0.023766751983212247, 'mean_word_length': 4.3962426010122675, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Search for Wandering Black Holes in the Milky Way with Gaia and DECaLS

We present a search for"hyper-compact"star clusters in the Milky Way using a combination of Gaia and the Dark Energy Camera Legacy Survey (DECaLS). Such putative clusters, with sizes of ~1 pc and containing 500-5000 stars, are expected to remain bound to intermediate-mass black holes (Mbh~10^3-10^5 M-sun) that may be accreted into the Milky Way halo within dwarf satellites. Using the semi-analytic model SatGen we find an expected ~100 wandering intermediate-mass black holes with if every infalling satellite hosts a black hole. We do not find any such clusters in our search. Our upper limits rule out 100% occupancy, but do not put stringent constraints on the occupation fraction. Of course, we need stronger constraints on the properties of the putative star clusters, including their assumed sizes as well as the fraction of stars that would be compact remnants.


WANDERING BLACK HOLES
Supermassive black holes are ubiquitous in the centers of massive galaxies today. They may play an important role in regulating star formation in galaxies (e.g., Silk & Rees 1998;McConnell & Ma 2013;Kormendy & Ho 2013). We do not know when or how supermassive black holes are formed. Upcoming gravitational wave experiments like Laser Interferometer Space Antenna (LISA; e.g., Amaro-Seoane et al. 2015) will potentially be sensitive to the mergers of the first "seed" black holes. However, with limited knowledge of both the rates and the mass functions of these black holes, the LISA observations cannot uniquely determine either. An important complementary clue will be provided by the study of relic intermediate-mass black holes, with the mass distribution and environment of these black holes today carrying some memory of when and how they were formed (e.g., Volonteri 2010;van Wassenhove et al. 2010;Ricarte & Natarajan 2018;Bellovary et al. 2019). The * NASA Hubble Fellow † NSF Astronomy & Astrophysics Postdoctoral Fellow black holes that are found outside of galaxy nuclei may be particularly sensitive to the seeding mechanism (see arguments in Greene et al. 2020). We expect to find black holes wandering in galaxy halos regardless of seeding mechanism because by hierarchical merging, galaxies are accreted by the Milky Way throughout its history, and some will be totally stripped apart from a stellar nucleus (e.g., Zinnecker et al. 1988). Some of these satellites, at least at the massive end, may bring intermediate-mass black holes into the galaxy with them (Volonteri & Perna 2005;Bellovary et al. 2010). The central black hole can retain a small cluster of bound stars, which we will refer to as a hypercompact star cluster (e.g., Merritt et al. 2009;Lena et al. 2020). A related mechanism for depositing off-nuclear black holes is gravitational recoil of black holes from the centers of merger remnants at early times (e.g., Volonteri et al. 2003a;O'Leary & Loeb 2009). Finally, if black holes are made through gravitational runaway at star cluster centers (e.g., Miller & Hamilton 2002;Portegies Zwart & McMillan 2002), then there may be a population of ∼ 1000 M black holes whose surrounding cluster has mostly dissolved, with only a hyper-compact cluster remaining (e.g., Fragione et al. 2018). A primary goal of this paper is to search for these star clusters that would, in turn, point to a population of black hole wanderers.
We do see evidence that star clusters come into our galaxy with their satellite. More specifically, nuclear star clusters are massive stellar clusters that are found at the centers of galaxies. It is thought that some fraction of the most massive known star clusters are the remnants of a stripping process that leaves nothing more than the nucleus (see, e.g., Pfeffer et al. 2016;Neumayer et al. 2020). The cluster M54 is the nucleus of the disrupting Sagittarius dwarf (Ibata et al. 1995). Another likely stripped nucleus is ω Cen, which shows evidence for multiple episodes of star formation (e.g., Norris et al. 1996;Pfeffer et al. 2016). We have not yet found definitive evidence of intermediate-mass black holes in either cluster (Ibata et al. 2009;Noyola et al. 2010;Lützgendorf et al. 2012;Baumgardt 2017;Baumgardt et al. 2019), nor for a similar cluster (G1), in Andromeda (Gebhardt et al. 2005). So far there are no definitive detections of black holes in Milky Way globular clusters (see discussion in Greene et al. 2020).
There is, however, evidence in galaxies more massive than the Milky Way for a population of wandering black holes. There are numerous dynamical detections of black holes in "ultra-compact" dwarfs (Seth et al. 2014;Ahn et al. 2018;Voggel et al. 2019). A number of the most massive ultra-compact dwarfs now have dynamically detected black holes at their centers, with the black hole accounting for ∼ 10% of the mass of the system in some cases. As outlined above, the most likely explanation for these objects is that they were formed in a more massive galaxy that was then stripped when falling into its current host halo (e.g., Pfeffer et al. 2014).
Perhaps the most compelling case for a wandering black hole is the "hyper-luminous" X-ray source HLX1 (Farrell et al. 2009). HLX1 is an accreting black hole, with a likely mass of M BH ≈ 10 4 M (e.g., Davis et al. 2011;Webb et al. 2012). It is found in a cluster of stars (Farrell et al. 2014) sitting a few kpc from a more massive galaxy ESO 243-49 at z = 0.022, which is very likely the remnant of a stripped dwarf galaxy that was accreted by the more massive system.
Our goal here is to search for lower-mass stellar clusters that would be bound because of the presence of an intermediate-mass black hole. O' Leary & Loeb (2012) and Lena et al. (2020) investigate a search for hypercompact star clusters using stellar colors. In this work, we instead focus almost exclusively on spatial clustering information. As described in §2, we expect any clusters deposited from accreted satellites to be found relatively near to the Galactic Center (< 50 kpc) and thus to be rather large on the sky and easily resolvable with existing ground-based imaging surveys. Therefore, in §3 and §4, we describe a joint search using the spatial resolution of Gaia and the depth of the Dark Energy Camera Legacy Survey (DECaLS; Dey et al. 2019). In §5, we summarize the limits we derive from our non-detections and in §6, we consider the future prospects of this work.

EXPECTED PROPERTIES OF THE WANDERERS
Our goal is to search for the hyper-compact stellar clusters that should accompany an intermediate-mass black hole wandering through our galaxy. In this section we summarize the relevant theoretical understanding of the size, mass, and stellar content of these hypothetical objects, so that we can hone our search strategy. It is worth emphasizing that we tune our search parameters to the specific case of "ex-situ" wanderers that formed in an external dwarf galaxy and were subsequently accreted by the Milky Way. In order for the galaxy to be fully stripped, we expect such clusters to live relatively close ( 50 kpc) to the Galactic Center based on modeling presented in §2.2.
Numerous successful searches in modern wide-field imaging surveys for globular clusters, dwarf galaxies, and stellar streams have been carried out in the Milky Way over the past two decades (e.g., Willman et al. 2005;Belokurov et al. 2006;Koposov et al. 2008;Bechtol et al. 2015;Koposov et al. 2015;Drlica-Wagner et al. 2015;Shipp et al. 2018;Torrealba et al. 2019). However, no search has been tuned to the compact sizes and low numbers of stars that we believe the hyper-compact star clusters may contain. Therefore, we thought it worthwhile to perform a customized search matched to the small scales of our target population.

Predicted size and mass
The most detailed calculations of the dynamical evolution of bound remnant clusters have been made by O'Leary & Loeb (2009), followed by confirming simulations in O' Leary & Loeb (2012). These authors focus on wanderers formed through very early mergers, in which the merged remnant is ejected from the proto-Milky Way center via gravitational wave recoil (e.g., Peres 1962;Campanelli et al. 2007), a model also explored by Volonteri & Perna (2005). Their modeling of the subsequent dynamical evolution of the cluster is very likely relevant to the final mass and size of the star clusters we consider as well, although the mechanism for forming the wanderers is different.
If we assume that the black hole is able to retain roughly its mass in stars within its sphere of influence, then the calculations of O' Leary & Loeb (2012) suggest that by the present day, the cluster will lose roughly 60-80% of that mass through relaxation (see arguments in Rashkov & Madau 2014). At the same time, the clusters will grow in size as t 1/3 , leading to present-day clusters of 0.5-1 pc in size, with mass ∼ 20% of the black hole mass.
We highlight two major potential caveats here. First, there is the possibility that the compact cluster is dominated not by visible stars, but rather by stellar-mass black holes. Such compact remnants preferentially reside at the cluster core and they would not be visible by electromagnetic means (e.g., Baumgardt et al. 2019;Gieles et al. 2021). O'Leary & Loeb (2012) acknowledge this possibility as well, but do not pursue it. We will take the same approach, but also note that extreme mass-ratio inspiral events, the detection of the merger of a stellar-mass black hole with an intermediate-mass black hole, would be one way to detect such clusters in the future with LISA (e.g., Gair et al. 2010).
The second caveat is the possibility that there is no stellar cusp around the black hole. Bahcall & Wolf (1976) calculate the stellar distribution around a black hole that ensues when the black hole is embedded in a stellar cluster with much higher mass than the black hole. However, low-mass galaxies can have very low stellar densities if they are not nucleated, and while the nucleation fraction of ∼ 10 9 M galaxies is near unity, that fraction drops substantially at lower galaxy stellar mass (e.g., Sánchez-Janssen et al. 2019;Neumayer et al. 2020). Therefore, it is possible that in practice, the wandering black holes lack even the hyper-compact clusters that we consider here. To partially compensate for this possibility, we consider the case that only nucleated galaxies host observable wandering black holes in §2.2.

Predicted number density and radial distribution
There have been a number of estimates for the number of wandering black holes that we might expect to find in a Milky Way-like galaxy (e.g., Volonteri et al. 2003a). We update these predictions in two ways. First, we utilize updated scaling relations that account for both the available data and upper limits (Greene et al. 2020). There is considerable scatter in these relations, which are of course extrapolations to the mass regime of interest here. Nevertheless, new dynamical measurements and constraining upper limits from the centers of nearby low-mass galaxies (Neumayer & Walcher 2012;Nguyen et al. 2018Nguyen et al. , 2019a motivate us to recalculate estimates of wanderer number density. Second, considerable work on high-resolution hydrodynamical simulations have yielded a new generation of semi-analytic models that can simultaneously capture the evolution of satellites in the tidal field of their host, and allow us to examine a large suite of Milky Waylike models in contrast with individual high-resolution simulations that have been used for this purpose in the past (e.g., Rashkov & Madau 2014). Specifically, we use the semi-analytic model SatGen (Jiang et al. 2020). This model builds upon Monte-Carlo dark matter halo merger trees and follows the evolution of satellites using tidal-evolution tracks calibrated against high-resolution idealized N-body simulations. The model also takes into account the response of dark matter halos to baryonic feedback, as formulated from zoom-in cosmological hydro-simulations, and the gravitational influence of a galactic disk on satellite evolution. In this way, it is possible to generate a large suite of Milky-Way-mass systems emulating those from high-resolution zoom-in cosmological simulations regarding satellite statistics, with a more complete sampling of the halo-to-halo variance and at numerical resolutions comparable to or higher than that of the simulations.
We use SatGen models with present-day halo masses distributed uniformly in the range M halo = 10 12 − 10 12.3 M . These models include an evolving disk component and are made to emulate simulations of bursty stellar feedback such as the FIRE and NIHAO simulations. Satellites with halo masses before infall that are M h > 1.3 × 10 6 M are considered.
From the model, we are particularly interested in satellites that dissolve, which is assumed to happen in SatGen when the sub-halo is stripped to a mass of 10 6 M . As shown by Jiang et al., a large fraction of the most massive satellites, particularly those accreted early, do pass close enough to the host halo center to be dissolved. Of course, it would be very interesting to follow these systems dynamically, including the impact of a black hole and possible stellar cluster on the subsequent stellar content and dynamical friction, but that is beyond the scope of the present work. For now, we simply count the number, infall mass, and radial distribution of these disrupted satellites ( Figure 1 & Figure  2).
Dissolved halos are still tracked as point masses by SatGen. We exclude halos that are within 10 pc of the galaxy center, assuming these effectively have merged with the central black hole. In practice such occurrences are very rare, although the satellite orbits are drawn from a cosmological distribution neglecting mass, so there is some chance that SatGen underestimates can- The cumulative number of black holes predicted assuming that the only channel to make wanderers is through accreted satellites. We estimate the central MBH of each satellite based either on the stellar mass (teal) or velocity dispersion (blue) at the time of accretion; the scatter includes both the range of accreted satellite masses in the SatGen models and scatter in the scaling relations. We show the 1σ range. As these are both extrapolations, they predict very different mass functions for the same population of disrupted satellites. We assume 100% occupation fraction in this panel, which is thus an upper limit on the number of wanderers formed via accretion of satellites. Roughly speaking, the clusters should contain ∼ 20% of the mass of the black hole nibalism if more massive satellites are preferentially on more radial orbits.
To convert the observed halo and stellar properties at infall into estimated black hole masses, we must decide what fraction of satellite halos will host black holes. The so-called "occupation fraction" of black holes in dwarf galaxies is not yet well measured. However, given the upper limit on a black hole in the nearby dwarf galaxy M33 (Gebhardt et al. 2001), it is likely that not every dwarf galaxy hosts a massive central black hole. We first calculate the predicted numbers of black holes under an assumption of 100% occupation, but additionally calculate the expected number of wanderers if we were to assume that only nucleated galaxies (galaxies with nuclear star clusters) host black holes. We adopt the nucleation fraction from Sánchez-Janssen et al. (2019), specifically their measurements in the Virgo Cluster. Virtually all galaxies with M * ≈ 10 9 M host nuclei, falling with mass such that galaxies with M * < 10 6 M will host black holes < 10% of the time.
We compute two possible black hole masses. First, we use the stellar mass at infall from the model combined with an M BH -M * relation from Greene et al. (2020). There are two major uncertainties here. First, we have assumed a stellar-to-halo mass relation to assign each halo a stellar mass (in this case the relation from Rodríguez-Puebla et al. 2017). The relationship between stellar and halo mass is notoriously unconstrained at these dwarf masses, with large degeneracy between the scatter and slope of the relation (e.g., Munshi et al. 2021), leading to large uncertainty in the stellar masses within the SatGen model. Second, the relationship between stellar mass and black hole mass has considerable intrinsic scatter and is dependent on galaxy morphology at higher mass (e.g., Reines & Volonteri 2015). We are extrapolating a high-scatter relation into an unknown regime, which adds considerable systematic uncertainty to these estimates.
As a second M BH estimate, we take the circular velocity at infall and, following Rashkov & Madau (2014), we calculate σ * = v circ /2.2 and then use the M BH -σ * relation (Greene et al. 2020). In this case, the circular velocity is securely predicted from the model, but the conversion between circular velocity and stellar velocity dispersion is not well known in the dwarf regime. As far as the scaling relations, the M BH − σ * relation is also somewhat morphology dependent (e.g., Greene et al. 2016), but the scatter is lower than the conversion based on total stellar mass. We will take the v max -based predictions presented below as the primary predictions 10 1 10 2 R (kpc) 10 0 Radial distributions of clusters that host black holes with MBH> 10 3 , 10 4 M where black hole mass is estimated following the MBH − σ * relation (blue-dotted and purple-solid respectively) or the MBH − M * relation (teal long-dash, all black holes with MBH> 10 3 M ). We also schematically indicate the typical radius of black holes made in recoil events from early mergers or from accreted satellites from Volonteri & Perna (2005). We highlight that wanderers deposited by disrupted globular clusters would be very close to the Galactic Center (Fragione et al. 2018), and we note that O'Leary & Loeb (2012) do not present a radial distribution, but do predict that black holes will extend to 100 kpc.
throughout the paper, because at least we will not be directly dependent on an assumed stellar-to-halo mass relation. At the end of this section, we will briefly discuss a third model for black hole mass in which there is no scaling with galaxy properties and all black holes are relatively massive at birth.
In calculating the expected number of wandering black holes, we include the halo-to-halo scatter by using the 85 Jiang et al. Milky Way-like models. We also vary the mapping between halo and black hole mass, by drawing 100 black holes per satellite from within the published scatter in the scaling relations. The resulting cumulative numbers of predicted wanderers are shown in Fig. 1 for the σ * and M * scalings respectively.
Roughly speaking, the stellar mass in the clusters will be ∼ 20% of M BH . Thus, we limit our attention to black holes > 10 3 M , where we might hope to detect the bound stellar cluster. Under an assumption of full occupation, we expect to find a few to 100 clusters with a hundreds to a few thousand stars, depending on both the black hole scaling relation adopted and on the oc-cupation fraction of black holes in the infalling satellite population.
We have an estimate of the satellite stellar masses from a stellar-to-halo mass relation (Rodríguez-Puebla et al. 2017). Assuming black holes with M BH > 10 3 M , in the case of the M BH -M * relation, nearly all hosts will have M * > 10 7 M . In contrast, in the σ * -based scaling, much lower-mass satellites host black holes. Specifically, we find that galaxies with stellar mass > 10 5 M are predicted to host black holes when scaling with σ * . Thus, we find many more possible black holes in the σ * -based scaling. Also, we note that even if we were to extend the lower-limit on halo masses considered by the model, we would not find any more black holes with M BH > 10 3 M .
Depending on the scaling relation, we not only expect different stellar masses, but also a different radial distribution, of satellites hosting black holes, as shown in Fig. 2. The most massive satellites, those that preferentially host black holes in the M BH -M * case, must be accreted early and travel close to the galaxy center to be stripped (Jiang et al. 2019). In contrast, the much wider range of halo and stellar mass that can host black holes when scaling with σ * also translates to a wider radial range, R 50 kpc. We will use the larger distance limit in §5 when we calculate upper limits on the number of wanderers in the Milky Way.
In Figure 1, we compare with existing similar predictions for the number of wandering black holes in a Milky Way-mass halo. Each paper makes slightly different assumptions about seeding mechanisms, and thus predict different mass spectra for the resulting wanderer population. Volonteri & Perna (2005) use a semi-analytic model as well, and directly consider multiple seeding mechanisms, to predict between one and ten black holes with M BH > 10 3 M . Bellovary et al. (2010) predict slightly higher numbers (5-15), in this case with a heavy seeding model. Rashkov & Madau (2014), like us, have no seeding model, but instead assign M BH based on properties of the halos, and predict comparable cumulative numbers as we do, using a similar M BH -σ * relation based on peak maximum halo velocity. Both Tremmel et al. (2018) and Ricarte et al. (2021) find ∼ 5 − 10 wandering black holes in Milky Way-mass halos using the Romulus simulation. Our predictions, and those from the literature, span two orders of magnitude in number because we are extrapolating assumptions about black hole scaling relations and occupation fractions into an unknown regime. These are the primary systematic uncertainties in our predictions, rather than detailed assumptions in the cosmological models that we use.
The predicted radial distributions in prior work are also similar to ours (such as Bellovary et al. 2010;Rashkov & Madau 2014), particularly those that model wandering black holes from disrupted satellites alone. The distribution is more extended than the wandering black holes expected from the dissolution of globular clusters should black holes form efficiently in their centers (e.g., Fragione et al. 2018), but is more centrally concentrated than models like Volonteri & Perna (2005) or O'Leary & Loeb (2012) that also include a component from early recoil events. In this work we do not attempt to search for more distant (spatially unresolved) clusters, which would require a different approach.
Finally, it is interesting to consider how our predictions would change if we adopted a heavy seeding prescription with a lower-mass limit of 10 4 − 10 5 M , as might be expected in some heavy seeding models (e.g., Bellovary et al. 2010;Inayoshi et al. 2020). If we assume that every satellite, regardless of mass, is seeded with a heavy seed, then the number of expected black holes with M BH > 10 3 M would be very high (∼ 1000) and would become very sensitive to our halo mass limit. Such a model is very easily ruled out, as we will show in the following sections. If, on the other hand, we assume a reasonable drop in occupation fraction with mass, then the number and mass distribution of black holes under a heavy seeding model would be very similar to Figure 1 (right), since the occupation fraction becomes a strong function of stellar mass, and this factor matters more than exactly how we assign black hole masses.

Predicted stellar content
We use the ArtPop software package (Danieli et al. 2018;Greco et al. 2021;J. Greco & S. Danieli, in preparation) to generate model clusters for calibration and completeness analyses. Using the MIST isochrones, we synthesize stellar populations of a given age and metallicity, and simulate realistic images of stellar systems based on their physical and structural parameters. In this work, we use ArtPop models to generate colormagnitude-position predictions in Gaia and DECaLS for expected clusters.
To build our model clusters, we assume a standard Kroupa (2001) initial mass function, with slope ψ ∝ m −α with α = 1.3 for stars below 0.5 M , and α = 2.3 at higher stellar mass. We set the lower and upper mass limits based on the MIST isochrones of a given age and metallicity. As described in §2.1, we do not explore the possibility that the present-day mass function is dominated by stellar-mass black holes. We assume a Plummer (1911) profile with a scale radius of 1 pc; the clusters are expected to grow to this large size with time (O'Leary & Loeb 2012). We assume a fixed size for our mock clusters. If in reality they are considerably larger (or smaller) by more than a factor of two, then our search is unlikely to find them.
We assume the clusters are 10 Gyr old, similar to their likely accretion time. To span the possible metallicity range of such clusters, we take [Fe/H]= −1.5, −0.5, which is measured for nuclear star clusters in host galaxies with M * = 10 8 − 10 9 M (Neumayer et al. 2020). In practice, the nuclei of such low-mass galaxies accreted at z > 1 are likely to have even lower metallicities, but we conservatively adopt these limits since the hotter and brighter stars in clusters at lower metallicity are easier to find at fixed stellar mass.
In Figure 3, we show the color-magnitude diagrams (bottom panel) and spatial distributions (middle panel) of example model clusters at the limits of our detection threshold ( §4.3), with stars detectable in Gaia and DE-CaLS in red and grey, respectively. The upper panel shows the artificial clusters injected into DECaLS images. It is worth noting that there is a large amount of stochasticity in the number of giant stars for a N = 1000 star cluster. While the median number of expected stars is five, it can range from zero to 15. We are not complete beyond the distance where we lose main-sequence stars ( §4.3) for these low-N models, and so we are not very sensitive to this stochasticity here. But, future searches at larger distances would be unless they are deep enough to reach below the main-sequence turn-off.

DATA
Having defined the parameters of the expected hypercompact clusters surrounding the putative intermediatemass black holes, we now turn to search for them. We first describe the data sets, and then in §4 we describe our search.

The Gaia catalog
Gaia is a European Space Agency astrometry mission. Here we use the Early Data Release Three of the Gaia mission (Gaia EDR3; Gaia Collaboration et al. 2016Collaboration et al. , 2018Collaboration et al. , 2020 to make an initial search for stellar overdensities on ∼ 15 scales by looking for stars with anomalously large numbers of neighbors on this scale. To do this, we utilize a custom-built catalog from S. Koposov 1 based on the EDR3 that records the number of neighbors that each source has with angular separation less  . These clusters should be comfortably detected by our Gaia search (see §3). Middle: The spatial distribution of the stars in these clusters as seen by Gaia (red) and DECaLS (grey). Bottom: Noiseless color-magnitude diagram for the clusters, colors as above. We should note that the number of giant stars for these low-mass clusters will vary quite significantly, from near zero up to fifteen for the N = 1000 star cluster, with a median of five. than x arcsec for 10 aperture choices of x. The aper-tures range from 0.25 to 128 arcsec, with each aperture increasing by a factor of two. It is convenient to work with the counts within "annuli" (i.e. number of neighbors between y and x arcsec), since each annulus can be modeled independently. We will refer to the number of neighbors between x/2 and x arcec as annuli counts A x . Our fundamental data set is therefore a vector of counts within these apertures for each EDR3 source.
3.1.1. Masking of known galaxies and star clusters There are two main types of contaminant that we mask before performing our cluster search. In general we mask a two degree radius around each satellite, except in the cases of the Large and Small Magellanic Clouds, along with Fornax, which we mask with a 10 degree radius (while the Large and Small Magellanic Clouds themselves do not fall within the DECaLS footprint, their outskirts turn out to be a major contaminant).
We also find that the cores of nearby galaxies can appear in Gaia as a set of point sources, and so we also mask the New General Catalog (NGC) galaxies 2 . The majority of these galaxies are given a 1-2 radius mask, while a small subset of galaxies with larger tabulated radii are masked over a full degree radius. In all, we mask roughly 1000 deg 2 of the 9000 deg 2 area covered by DECaLS.

DECaLS
The DECam Legacy Survey (DECaLS) has been carried out with the Dark Energy Camera (DECam; Flaugher et al. 2015;Dey et al. 2019) at the Mayall 4m telescope at the Cerro Tololo Inter-American Observatory. DECaLS reaches 5σ point-source depths of grz = 23.95, 23.54, 22.50 AB mag over 9000 deg 2 . We utilize gr photometry from Data Release 9 of DECaLS 3 in our search.

SEARCH WITH GAIA+DECALS
As described in §2, we seek hyper-compact star clusters comprising 500-5000 stars within ∼ 1 pc. These clusters likely fall within 50 kpc of the Galactic Center (see Fig. 2). In the proof of concept search that we present here, we have decided to combine the spatial resolution of Gaia with the depth of DECaLS to search efficiently for possible clusters. First, we use the Gaia EDR3 catalog to search for clusters of N > 4 stars within 15 , as expected for these hyper-compact star clusters. The benefits of Gaia include the high spatial resolution and the potential to filter on proper motions. However, with a depth of G ∼ 21.5 mag (Gaia Collaboration et al. 2020), we will detect only a few of the stars with Gaia (as we see in Fig. 3).
To illustrate this point, in Figure 4, we show the expected number of stars that can be observed with Gaia and DECaLS in a cluster of 500-5000 stars as a function of distance. We use 100 different realizations of the cluster to include stochasticity in how many stars populate the red giant branch for the lower-mass clusters. For N = 1000, we expect a large range in the number of giant stars, between zero and twelve, with a median of five. However, over the distance range where we are sensitive with this Gaia search ( §4.3), we do reach below the main-sequence turn-off, meaning that our overall N star,Gaia is dominated by main-sequence stars.
The grey regions show the numbers for the highermetallicity clusters, showing that the range in metallicity leads to only a small difference in the number of stars, so throughout we will focus on models with [Fe/H]= −1.5, which seems more likely. We see that overdensities of stars identified with Gaia should be accompanied by an increase of three to five times in the number of stars at detected at DECaLS depth, if we are identifying star clusters. This will be true even with some mild crowding, as discussed in §4.2.1.
In this section, we first present the model that we use to identify possible clusters as outliers in spatial clustering ( §4.1). We then calibrate the methods using artificial clusters ( §4.2), argue that crowding is unlikely to greatly impact our search ( §4.2.1), and determine our completeness ( §4.3). Finally, with a list of cluster candidates in hand from Gaia, we use DECaLS to invalidate most candidates based their total number of stars at DECaLS depths ( §4.4) along with their colors ( §4.5).

Negative Binomial Model
Our method relies on modeling the distribution of neighbor counts around each star in annuli (as described in 3.1) and looking for outliers, or low probability points, in this distribution.
For each star, we model the vector A, which is the neighbor counts in angular annuli on the sky over all scales up to 128 . We model each annulus A x in a small region of sky as an independent probability distribution, with no correlations between annuli. Then we have that

Number
Starfield 500, 30kpc 1000, 50kpc 5000, 60kpc Figure 5. Log of the probability that the model stars (colored lines) are drawn from the parent distribution that describes the probability of the annuli Ax for a random field (black). The probabilities are calculated assuming that the counts within annuli are randomly distributed according to a Negative Binomial distribution.
where φ x are the model parameters. If all of the counts were random, we could model each probability distribution as a Poisson distribution. In practice, however, the dispersion may well be larger than the mean due to density variations over the fields of view that we consider.
In order to account for this possibility, we specify P (A x ) with the Negative Binomial distribution, which is a generalization of a Poisson distribution, and can be thought as a Poisson distribution whose mean is sampled from a gamma distribution. The Negative Binomial is a discrete probability distribution describing the number of 'unsuccessful' trials k that occur in some series of repeated, independent trials, each with success probability p, before r successful trials occur. This distribution can be written as This distribution approaches the Poisson distribution if one sends p → 1 and r → ∞ while keeping the mean of the distribution [λ = (1 − p)r/p] constant. Unlike in the Poisson distribution, the variance [σ 2 = (1 − p)r/p 2 ] and mean are two independent quantities. We determine the parameters of the Negative Binomial φ = r x , p x by modeling the A x of stars in neighboring regions set to be a few times the expected scale of the hyper-compact clusters themselves. We divide the sky into healpix pixels with NSIDE= 256 (∼ 2 ) pixels, and fit the Negative Binomial coefficients to stars in the surrounding nearest eight healpix pixels. We then use that parameter vectorφ, comprising each fitted r x , p x for each annulus, to estimate the likelihood of each star P (A|φ). We find p to range from 0.7-0.9 at large radii, with significant improvements to χ 2 over a Poisson fit. In the smallest annuli < 8 , p approaches 1, and we therefore adopt the Poisson distribution fits at these smallest annuli.
The final probability is calculated as the product of the probabilities from the fit to each annulus (Equation 2 above). Stars that live in anomalously dense regions of sky on ∼ 15 arcsec scales will have very low probability in this model, which assumes that stars are distributed randomly. In the next section we present mock tests that we use to pick an outlier threshold.

Selecting outliers
To hone our outlier selection and test our completeness, we use the ArtPop models described in §2.3 to inject artificial clusters into our data, and calculate their Negative Binomial probabilities, as a function of number of stars in the cluster, distance to the cluster, and the background stellar density in each field.
We have generated 100 realizations for each of three cluster masses (N star = 500, 1000, 5000 M ) and two metallicities ([Fe/H]= −1.5, −0.5). These 100 models are each inserted into 100 random locations that span a range of background stellar density, for a total of 10 4 models for each N star . For book-keeping purposes we track only the star closest to the center, creating annular counts that include both artificial stars in the cluster and real stars in the random location we have chosen. We then calculate the probability of each star in the Negative Binomial model. While these 100 models capture the variation in stars on the main sequence well, they will not capture the full stochasticity on the giant branch. Therefore we ran a second set of 1000 realizations for the N = 500, 1000 low-metallicity models to verify that our completeness and thresholds do not change.
The distributions of probabilities from the Negative Binomial fits are shown in Figure 5 for the lowermetallicity [Fe/H]= −1.5 model. For display purposes, we select models at a distance and number of stars that is roughly at our 80% completeness limit, as shown in Figure 6. We then select a probability threshold for candidate clusters that retains a high fraction of the model cluster stars without swamping us with "normal" background stars. There are very few stars with probabilities log P < −10, but we can still maintain high completeness in the mock clusters ( §4.3). Generally, with this cut we select a fraction 0.01% of the stars, for a total of 22,200 candidates.
In principle, we also have access to proper motion data for the Gaia stars. We would expect the hyper-compact clusters to have velocities and positions consistent with halo objects, and we would expect the stars to have proper motions consistent with each other, as the inter-nal motions of the cluster stars should be smaller than their bulk motions. However, in practice we are finding typically only 3-4 stars that may be associated with each other ( §4.5), and some of these are too faint for reliable proper motions. In practice we therefore do not use the proper motion measurements in our search.

Possible impact of crowding
Here we investigate what fraction of stars in the core of our clusters may be lost in the DECaLS catalog due to crowding (Fig. 3). While the clusters we are seeking will be a few arcsec across on the sky, as Figure 3 shows, we could expect some blending to affect the number of stars in the DECaLS catalogs. While this blending will not impact our completeness corrections, that are based entirely on Gaia, they could impact our ability to detect the larger population of cluster stars that we expect at DECaLS depths.
We expect there to be 0.07 stars arcsec −2 within the half-light radius of a 500-star cluster at our detection limit of ∼ 30 kpc, 0.17 stars arcsec −2 in a 1000-star cluster at our limit of ∼ 40 kpc, and 0.9 stars arcsec −2 in a 5000-star cluster at 50 kpc, beyond which we do not expect many clusters. We investigate how the DE-CaLS catalog behaves using the star clusters Reticulum and Whiting 1. Denser and more luminous clusters are modeled differently within the DECaLS catalog, which uses Gaia data to identify stars and does photometry only on these stars 4 .
To get a handle on crowding using these clusters, we make two assumptions. First, we assume that all the objects at the centers of these two clusters are actually stars (probably nearly true). Second, we consider each object of any type to represent a single star. This is not always true. For instance, a single exponential object can be a blend of two or three stars. However, it will at least give us an idea of the magnitude of crowding. Whiting 1 has 0.09 stars arcsec −2 , with 70% of those identified as point sources. Reticulum is denser, with 1.7 stars arcsec −2 and 50% of those identified as stars. Thus, at the crowding levels that we expect for the majority of our cases, it is safe to assume that ∼ 30% of stars may be lost to crowding. Note that things will be worse in the center of a 5000-star cluster, but there would be many more stars outside the core that we would find. This is of course approximate, but it demonstrates that at the level of crowding that we expect, it is still safe to use the DECaLS photometry. We also note that we do not expect any saturated stars.  Figure 6. The fraction of clusters that we recover as outliers from our Negative Binomial fit for three cluster masses (500 stars in blue, 1000 stars green, 5000 stars red) as a function of distance. We define our survey volume by the distance where we reach 80% completeness (dashed vertical blue for 500 stars, dashed vertical green for 1000 stars). Based on the SatGen models, we expect most hyper-compact clusters would lie within 50 kpc of the Galactic Center, excluding the grey shaded region.

Completeness
We use the ArtPop simulations to quantify our completeness. For each cluster mass, and each distance, we first verify that there is no difference in completeness as a function of background stellar density in each field. This independence is not surprising, given that the search for very compact objects is limited by Poisson noise in the number of stars at the cluster scale rather than the Poisson noise in the background. We then measure the fraction of clusters that are identified by our outlier threshold based on the Negative Binomial fits, as a function of distance in Figure 6. Our probed volume is defined within the distance where we fall to 80% completeness (30, 40 kpc for 500, 1000-star clusters, and beyond 50 kpc for 5000-star clusters). We will use these numbers to quantify the limits on the number of such clusters that may still be lurking in the Milky Way in §5.

Eliminating candidates with DECaLS matching
As shown in Figure 4, we expect to detect many tens of stars from a cluster at our detection limit in DECaLS, where we might only expect four to six stars in Gaia. Therefore, a straightforward way to winnow down candidates identified with Gaia is simply to ask how many stars are detected at that position in DECaLS. From the DECaLS photometry, we count the number of stars and extract color-magnitude distributions, both of which can be used to determine whether or not we have found hints of a real cluster.
By cross-matching with DECaLS, we are able to apply the following additional cuts. We only consider candidates that have at least N star,DECaLS = 9 stars in the DECaLS catalog within 15 arcsec of the star in question. The aperture is chosen to be roughly twice the effective radius of a cluster around our completeness limit. We also apply a color criterion that the stars must fall in the range 0.2 < g − r < 1.5, inspired by the color range of the color-magnitude diagram of an old and metal-poor stellar population (see Figure 3). Selecting nine stars as a limit is conservatively lower than the number we expect within < 15 arcsec (25 ± 5 at our detection limit for a 500-star cluster, Fig. 4). Thus we allow for some loss of stars from confusion. However, with nine stars we still have sufficient numbers to crudely fit a color-magnitude relation ( §4.5).
As an additional sanity check on the impact of crowding, we also include in our "star" count objects that DECaLS has classified as "compact" exponential sources (REX). We only consider REX sources with sizes smaller than the PSF, but in this way we can crudely evaluate whether we are missing stars that have been classified as extended due to blending. In practice, we find that adding these makes little difference to our final list of targets, since most REX-dominated candidates tend to be distant galaxy clusters that are eliminated by our color cut.
We find 176 stars from 86 candidate clusters with at least nine DECaLS stars falling in a reasonable color range, once we remove some contamination from inadequately masked nearby galaxies.
We expect ∼ 25 ± 5 stars in a 15 aperture for the 500-star clusters at our detection limit (Fig. 4). There is only one candidate cluster in our outlier sample with N star,DECaLS > 12, and 6 with N star,DECaLS > 11. This number is already ∼ 2σ lower than the low end of the expected range at our distance detection limit, assuming Poisson errors. However, in order to determine whether we have identified any candidates, it is useful to also examine the colors and magnitudes of the candidates. With color and magnitude we can isolate associated stars from foreground/background objects.

Eliminating candidates with a color-magnitude fit
For all cluster candidates, we next determine whether the DECaLS gr color-magnitude distribution is consistent with an old coeval stellar population. We limit our attention to those candidates that have at least N star,DECaLS = 9 stars with 0.2 < g − r < 1.5 within 15 . An example of the isochrone fit to the candidate cluster shown at left. The red points are drawn from a 30 radius around the candidate, while the grey-scale represents our best-fit isochrone. Even with this wider aperture, we still only find a maximum of six stars that may be associated with the same isochrone.
To fit the color-magnitude diagrams (CMDs), we include two components, a background model built from the data, and a simple stellar population model at a range of distances. The model of the background is constructed from all of the stars in 10 degree patches. The possible additional cluster component is made from the MESA Isochrones and Stellar Tracks isochrones (MIST; Dotter 2016; Choi et al. 2016) built from the Modules for Experiments in Stellar Astrophysics (MESA Paxton et al. 2011) models sampled with a Kroupa IMF. For simplicity, and because we do not expect young stars in our hyper-compact star clusters, we fix the age to 10 Gyr and the metallicity to [Fe/H]= −1.5. There are two free parameters in our fit, the relative fraction of background and cluster stars, and the distance for the cluster component, allowed to range from 10 to 90 kpc.
We build a two-dimensional histogram of the theoretical color-magnitude diagram thus constructed, including a scatter in the magnitude and color derived from the DECaLS data. We then select the distance modulus and background weight that maximizes the likelihood (log L = Σ i logCMD i ), where CMD is the normalized probability density at the position of each of i stars in our candidate cluster. An example fit is shown in Figure  7.
We test this method using our mock clusters generated by ArtPop, described above in §2.3. Specifically, we embed the model clusters in random DECaLS fields, add photometric noise, and then run our fit on this collection of stars. Once main sequence turnoff stars are too faint to be detected by DECaLS, we can no longer estimate a reliable distance (roughly > 40 kpc), but we can achieve reasonable fits when the number of stars in the cluster is N star,DECaLS > 9.
Our most likely CMD fits yield a refined N star,DECaLS for each candidate cluster. Although we have identified overdensities in space, the colors of the stars suggests that the majority are foreground/background stars. According to our fits, the largest number of stars that may be associated with a common isochrone is N star,DECaLS = 4 − 5 stars (two candidates), or N star,DECaLS = 3 stars, (14), while all the other potential cluster stars are drawn from the background. Since the low end of our predicted number of stars is 25 ± 5 at our detection limit, we conclude that with high significance, we have not detected a candidate hyper-compact star cluster as defined in §2.
As an example of an apparent overdensity, we highlight one candidate cluster from the list of 86 ( Figure  7). Within 15 , this candidate has three stars associated with the same isochrone, which grows to 5-6 when we open the aperture to 30 . This candidate comprises a clear overdensity spatially. As a sanity check, we fit Poisson distributions to randomly selected 15 apertures within 3 of the candidate, and find that it is a 3 σ outlier in density. However, as judged from the isochrone fit, only a few of those stars may be associated with each other. They do have consistent proper motions, but even if they are physically associated, there are not enough of them to represent the clusters that we are looking for.

LIMITS ON WANDERING BLACK HOLES IN THE MILKY WAY
We searched ∼ 8000 deg 2 of sky for hyper-compact star clusters of 500-5000 stars within 50 kpc of the Galactic Center, that would be the signature of wandering black holes accreted with infalling satellites. We did not find any plausible candidates, so in this section we translate our non-detection into upper limits on this population.
To calculate the upper limits, we assume that we are sensitive within the volume defined by our 80% completeness, as indicated in Figure 6. We further assume that the black holes lie within 50 kpc, as calculated using the M BH − σ * relation to populate halos with black holes (see Figure 2).
The resulting limits are shown in Figure 8 and Table 1. The predicted number of black holes (or corresponding star clusters on the top axis) under the assumption of 100% occupation fraction are shown by the dotted line, while taking the nuclear cluster occupation fraction as the black hole occupation fraction is shown in solid. In what follows we discuss what we can and cannot conclude from these limits.

Limits on the occupation fraction
In this paper we have searched for and place limits on the number density of star clusters with < 5000 stars. We are interested in the corresponding number of black holes. We translate between the cluster mass limit and M BH by simply assuming that the cluster mass is 20% of M BH .
The number of black holes that we expect is basically the product of the halo mass function, the occupation fraction, and a convolution with the halo to black hole scaling relation. If we assume that σ * can be related to the maximum halo circular velocity as σ * = v max /2.2, then we can take the black hole scaling relation with v max to have the form: log M BH = α + βlog v max . For intrinsic scatter in this relation of σ int , we then have (see, e.g., Marconi et al. 2004;Gallo & Sesana 2019): We then write the number of black holes in terms of the halo maximum velocity function, the occupation fraction λ occ (v max ), and P (log M BH |v max ): where the halo velocity function is given by N h (v max ) and the occupation fraction is λ occ (v max ). The halo mass function in our predictions is derived from the SatGen simulation. The M BH scaling relation is extrapolated from the observed M BH − σ * relation, assuming that the maximum halo circular velocity can be converted directly to a stellar velocity dispersion with a constant scale factor. Both the extrapolation and the conversion are uncertain. As described in §2.2, we do not consider the M * −M BH scaling relation here, given the additional uncertainty added by the unknown stellar-tohalo mass relation.
Currently, the occupation fraction is not constrained in the regime of interest to us, below galaxy stellar masses of M * ≈ 10 9 M , or σ * ≈ 20 km s −1 . There are observational constraints on occupation fraction for more massive dwarf galaxies with stellar masses of M * 10 9 M (see summary in Greene et al. 2020), where the occupation fraction is consistent with at least 50% of galaxies hosting a central black hole, based both on stellar dynamical results (e.g., Nguyen et al. 2018Nguyen et al. , 2019b and X-ray studies (e.g., Miller et al. 2015;She et al. 2017). At stellar masses M * < 10 9 M , there are precious few constraints on the occupation fraction, although some AGN have been found in galaxies of this mass (e.g., Baldassare et al. 2019;Reines et al. 2019).
Our measured limits are inconsistent with 100% occupation for sub-halos of the Milky Way. They are well below predictions for the number of black holes we would expect if every satellite contained a black hole with mass as predicted by the M BH − σ * scaling. As mentioned in §2.2, a heavy-seeding model in which all halos down to ∼ 10 6 M are seeded is ruled out even more conclusively. To take a concrete example of a non-unity occupation fraction, we adopt the nucleation fraction (the fraction of galaxies containing nuclear star clusters) as a proxy for the black hole occupation fraction, shown as a solid line in Fig. 8). All black holes discovered dynamically in galaxies with M * ≈ 10 9 − 10 10 M are found in nuclear star clusters (e.g., Seth et al. 2008;Nguyen et al. 2018Nguyen et al. , 2019b, potentially suggesting a relationship between the two (see details in Neumayer et al. 2020). Nucleation fractions are near unity at ∼ 10 9.5 M , and then fall to < 10% by M * ∼ 10 6 M (e.g., Sánchez-Janssen et al. 2019;Neumayer et al. 2020).
Our observed limits are consistent with an occupation fraction that falls with mass in a similar way to the nucleation fraction. In fact, the details of how we seed the satellites become relatively unimportant in this case, because the number of black holes is set by the small number of available dissolved satellites with initial M * > 10 7 M , where most of the black holes would 5 × 10 2 10 3 24 10 3 5 × 10 3 11 5 × 10 3 10 4 5 10 4 5 × 10 4 3 Note-Limits on number of hypercompact star clusters, and corresponding limits on the number of wandering massive black holes inferred in this work. Theoretical expectations from SatGen, assuming either that every satellite hosts a 10 3 −10 5 M black hole are shown in dotted, or in solid the perhaps more realistical case that the fraction of galaxies hosting black holes is the same as the fraction of galaxies that host nuclear star clusters. Black hole masses are estimated from the halo maximum velocity (blue). Our limits assume that all clusters would be found within 50 kpc of the Galactic Center.
be found. If we improve the probed area by a factor of a few (e.g., with the Rubin Observatory; Ivezić et al. 2019), we might be able to constrain the occupation fraction mass dependence further. In the meantime, we still expect that there must be a wandering black hole population in Milky-Way like galaxies, and we explore other ways to search for them in §5.3.

Important limitations to our assumptions
Dynamical modeling that resolves the scale of the sphere of influence of the black hole is needed to be secure in the properties of the hyper-compact star clusters. There are at least three effects that we have ignored or oversimplified that could strongly impact our conclusions. First, we assume that once a satellite reaches a stellar mass of 10 6 M in SatGen, it will fully dissolve. This assumption should be tested in the context of a central intermediate-mass black hole.
Second, we have ignored the possibility that hypercompact star clusters may be composed mostly of compact remnants. Given that neutron stars and black holes will settle to the cluster center through masssegregation, there is some chance that these clusters exist but have much higher mass-to-light ratios than assumed here. In fact, a mass fraction in stellar-mass remnants of ∼ 20%, which sounds high compared with some models (e.g., Zocchi et al. 2019;Baumgardt et al. 2019), has recently been suggested to explain the low-density of Palomar 5 and potentially other fluffy globular clusters (Gieles et al. 2021 Third, we have assumed that all of the clusters have the same size of 1 pc, which is clearly an oversimplification and is also related to the detailed dynamical evolution of the cluster. If the clusters are significantly larger than the ∼ 1 pc that we assume here, then their clustered signal would drop, as would our sensitivity. As an example, if the N star = 1000 cluster had a size of 5 pc, then our completeness would drop roughly by a factor of two compared to what is shown in Figure 6.

Other ways to search for wandering black holes
While we have been focused on low-mass hypercompact clusters, one promising avenue for continued study is certainly to identify all of the Milky Way star clusters that may be stripped nuclei, using both stellar population and orbital information (e.g., Massari et al. 2019;Pfeffer et al. 2021). Future observations with extremely large telescopes will provide far more stringent limits on their possible central black hole content.
Another angle is to think about the other channels that may generate hyper-compact clusters surrounding black holes at larger galactocentric radius, including early ejection from the gravitational slingshot (e.g., Volonteri et al. 2003b;O'Leary & Loeb 2009;Lena et al. 2020). It is therefore interesting to consider searches that will reach larger galactocentric radii in the Milky Way, as well as complementary techniques for searching for extragalactic analogs.
Euclid and the Vera Rubin Observatories Legacy Survey of Space and Time (LSST; Ivezić et al. 2019) will provide a more sensitive search of these potentially more distant hyper-compact star clusters. Here we consider what these clusters would look like in next-generation imaging surveys. From the ground, the LSST will produce single-epoch images with 5 σ depths of r = 24.7 mag. At this depth, we could expect to uncover 500-star hyper-compact clusters out to ∼ 70 kpc, where we still expect ∼ 25 ± 5 stars to be detected within a 3 radius, although crowding will likely be significant. A 1000-star cluster will have r = 17 − 18 mag within ∼ 2 at 70-90 kpc, making them easily detectable in principle. However, these clusters will be very crowded at ground-based spatial resolution, and so techniques using colors may be required for these more massive clusters beyond ∼ 50 kpc. Such clusters may be selected as outliers in the color-magnitude diagram, as they will be far more luminous than expected for their color (e.g., Lena et al. 2020). One might imagine combining such a ground-based search with higher-resolution data from a mission like Euclid (Racca et al. 2016) to see if the putative clusters are resolved into a couple of stars at higher resolution.
Finally, an additional prospect is opened by the Roman Space Telescope (Spergel et al. 2015) to search for the ∼ 5000 star hyper-compact clusters in halos of external galaxies. While these clusters would be point sources, they would be anomalously bright for their color. Since such a cluster would have an integrated magnitude of z ≈ 22 − 23 mag at 10 Mpc, one could imagine performing a search leveraging the proposed Roman Infrared Nearby Galaxy Survey (RINGS) (Williams 2015) to search for hyper-compact star clusters at the very low-luminosity end of the star cluster luminosity function, using the high spatial resolution to remove background galaxies and the colors to distinguish from halo stars, as extragalactic globular clusters are currently found. However, ultimately identifying these as real cluster candidates will be very challenging, as will certainly require spectroscopy.
Another tool, and the most robust one, is to search for the dynamical signatures of black holes in the motions of stars. At somewhat higher mass, stripped remnant clusters surrounding massive black holes have been found from dynamical signatures in ultra-compact dwarf galaxies (Seth et al. 2014;Ahn et al. 2017Ahn et al. , 2018. Closer to home in Andromeda, lower-mass stellar clusters, likely to be ultra-compact dwarfs, also show some signs of nuclear black holes (Gebhardt et al. 2005) al-though these remain very challenging measurements. In the era of extremely large telescopes, we can hope to detect 10 4 M black holes dynamically out to ∼ 5 − 10 Mpc, which will reach a much larger number of massive globular cluster and UCDs .
There are promising search avenues via accretion signatures as well. One is through tidal capture of individual stars in the cluster (e.g., MacLeod et al. 2016), which may be what we are observing in the hyper-luminous Xray source HLX1 (e.g., Farrell et al. 2009;Lasota et al. 2011;Soria et al. 2017). Tidal disruption events should probe wandering black hole populations, and there are some potential tidal events in off-nuclear stellar systems (e.g., Lin et al. 2018).
It is possible that wandering black holes in more massive halos (e.g., the Virgo cluster) may encounter sufficient gas to be detected via accretion of inter-cluster medium (e.g., Guo et al. 2020). This detection channel is highly unlikely in Milky Way-mass systems. It may, however, be possible to detect very low Eddington-ratio accretion onto black holes in the centers of normal extragalactic globular clusters in the radio, which is relatively boosted at low accretion rates (e.g., Maccarone et al. 2005). On the other hand, searches with deep radio surveys have not been successful to date, either in the Milky Way (Tremou et al. 2018) or external galaxies (Wrobel et al. 2016). Next-generation radio facilities will be very powerful to further the reach of these searches (Wrobel et al. 2019).
Finally, wandering black holes may eventually be detected as extreme mass-ratio inspiral events through the merger of a stellar mass black hole with the central intermediate-mass black hole (e.g., Gair et al. 2010;Amaro-Seoane et al. 2015;Eracleous et al. 2019). Using the tools developed here, it would be useful to calculate the number densities of extragalactic wandering black holes to estimate a detection rate for space-based gravitational wave experiments.

SUMMARY
We have presented a search for hyper-compact star clusters that would be the hosts of wandering black holes. If such black holes are deposited by infalling satellites, we expect them to be ∼ 1 pc in size, have 500-5000 stars, and fall within ∼ 50 kpc of the Galactic Center. We use Gaia+DECaLS data over ∼ 8000 deg 2 . Using a Negative Binomial model to describe the distributions of stellar counts around each target star in Gaia EDR3, we identify large outliers in count space. Real stellar clusters with normal mass functions would have threefive times as many stars in DECaLS, allowing us to efficiently eliminate candidates through a cross-match with DECaLS. We do not find any hyper-compact star clusters within ∼ 30 − 50 kpc (for clusters of 500-5000 M ). We translate these limits into upper limits on the number of intermediate-mass black holes wandering in the inner Milky Way halo.
We also use modern semi-analytic models to bracket the number of 10 3 − 10 5 M black holes that might wander in the Milky Way halo using the SatGen code (Jiang et al. 2019). SatGen calculates the number of expected dissolved satellites as a function of galactocentric radius and halo/stellar mass at infall. We then extrapolate black hole-galaxy scaling relations to predict the range of black hole masses that might inhabit these halos. Based on our calculations, we expect as many as 100 black holes with M BH > 10 3 M based on scaling from the M BH − σ * relation, should every satellite carry an intermediate-mass black hole. In the context of this model, we can rule out that all satellite galaxies host central black holes from our non-detection. Our measurements are consistent with models in which most black hole hosts have stellar masses > 10 7 M .
In the near future, wide-area imaging surveys from the ground and space will open up new discovery space for hyper-compact star clusters in the Milky Way and nearby galaxies. Extremely large telescopes will enable dynamical detection of 10 3 M black holes in the Milky Way and 10 4 M black holes within the Local Volume. The Rubin Observatory LSST will increase the rate of detection of tidal disruption events, some of which may be in off-nuclear systems. Eventually gravitational wave detectors in space will be sensitive to the mergers of stellar-mass and intermediate-mass black holes in these clusters (e.g., Gair et al. 2010;Gallo & Sesana 2019). Thus, if there are black holes wandering in galaxy halos, we will begin to uncover them relatively soon.