Hypercompact stellar clusters: morphological renditions and spectro-photometric models

Numerical relativity predicts that the coalescence of a black hole-binary causes the newly formed black hole to recoil, and evidence for such recoils has been found in the gravitational waves observed during the merger of stellar-mass black holes. Recoiling (super)massive black holes are expected to reside in hypercompact stellar clusters (HCSCs). Simulations of galaxy assembly predict that hundreds of HCSCs should be present in the halo of a Milky Way-type galaxy, and a fraction of those around the Milky Way should have magnitudes within the sensitivity limit of existing surveys. However, recoiling black holes and their HCSCs are still waiting to be securely identified. With the goal of enabling searches through recent and forthcoming databases, we improve over existing literature to produce realistic renditions of HCSCs bound to black holes with a mass of 10$^{5}$ M$_{\odot}$. Including the effects of a population of blue stragglers, we simulate their appearance in Pan-STARRS and in forthcoming $Euclid$ images. We also derive broad-band spectra and the corresponding multi-wavelength colours, finding that the great majority of the simulated HCSCs fall on the colour-colour loci defined by stars and galaxies, with their spectra resembling those of giant K-type stars. We discuss the clusters properties, search strategies, and possible interlopers.


INTRODUCTION
General relativity predicts that a burst of gravitational waves (GW) is emitted during the coalescence of two compact objects (e.g. Fitchett 1983;Redmount & Rees 1989;Wiseman 1992). The first direct observation of such an event took place on 2015 September 14, marking the advent of a new era. Abbott et al. (2016) inferred that the emission originated from the merger of two black holes with masses 29 +4 −4 and 36 +5 −4 M . Asymmetries in the merging objects (i.e. different masses and spins) are expected to produce asymmetries also in the GW emission, leading to a net flux of linear momentum. To conserve linear momentum, the merging binary and the resulting object recoil accordingly. They get a "kick". The amplitude of such kick is largest for black hole (BH) binaries, and it must be computed numerically. When this became E-mail: d.lena@sron.nl possible, the expected amplitudes of order 10 2 km s −1 were confirmed (Pretorius 2005;Baker et al. 2006;Campanelli et al. 2006;Pretorius 2006). Moreover, much larger values, as high as 5000 km s −1 for special configurations, were also obtained (e.g. Campanelli et al. 2007;Tichy & Marronetti 2007;Brügmann et al. 2008;Rezzolla et al. 2008;Lousto & Zlochower 2011, 2013. Via the modelling of GW waveforms, recoil velocities have been estimated for GW170104 (Healy et al. 2018), and for GW150914 , where the inferred kick attains −1500 km s −1 . ate kick and is not ejected from its host galaxy, is the one expected to be the most common: the recoiling BH would oscillate about the galaxy center loosing energy via dynamical friction at each passage through the core of an early-type galaxy (e.g. Gualandris & Merritt 2008), or through the disk and bulge of a late-type galaxy (e.g. Kornreich & Lovelace 2008;Blecha & Loeb 2008); bursts of accretion would follow the passage through a gas-rich disk (e.g. Blecha & Loeb 2008). The second scenario, where recoiling velocities are high enough to remove a SMBH from its host galaxy, is much more unlikely than the first: attaining recoiling velocities in excess of a few hundreds km s −1 requires special configurations for the BH-binary (e.g. Lousto & Zlochower 2011;Lousto et al. 2012). For instance, gas accretion onto the merging SMBHs could align their individual spins with the binary orbital angular momentum, heavily hampering the recoil velocity, and producing kicks which are, for the most part, below 100 km s −1 (e.g. Bogdanović et al. 2007;Dotti et al. 2010).
Depending on the environment where it is born, the recoiling BH carries along a mixture of gas and stars, with the amount of matter bound to the recoiling BH being inversely proportional to the kick velocity, V k . When the recoiling BH originates in a gas-rich environment (the "wet" merger scenario), for example from a binary embedded in a gaseous disk, then, assuming a disk with a mass much smaller than the binary, the recoiling BH is expected to carry with it a punctured disk with outer radius r ∝ V −2 k , and mass M disk ∝ V −2.8 k . This gas will be accreted within a few million years (e.g. Loeb 2007;Bonning et al. 2007).
When gas is not present (the "dry" merger scenario), then the recoiling SMBH will still carry with it a retinue of stars: those located within a distance r k ≡ GM • /V 2 k from the SMBH will remain bound after the kick, and the predicted stellar mass of the cluster is M ∝ V −2(3−γ) k , with γ the slope of the stellar density distribution before the kick (Komossa & Merritt 2008;Merritt et al. 2009;O'Leary & Loeb 2009). The effective radius for these clusters is predicted to depend on a number of parameters: the central velocity dispersion of the host galaxy, the stellar density distribution prior to the kick, and the dynamical status of the nucleus (collisional or collisionless, Merritt et al. 2009); while the largest clusters could extend as much as a few tens of parsecs (similarly to globular clusters and ultra-compact dwarf galaxies), the great majority is predicted to have sizes below 1 pc (hence the appellative "hypercompact") and velocity dispersion in excess of a few tens of km s −1 (Merritt et al. 2009;O'Leary & Loeb 2009), much larger than the typical velocity dispersion of globular clusters (approximately 10 km s −1 , e.g. Pryor & Meylan 1993).
The observed velocity dispersion of an hypercompact stellar cluster (HCSC) is of primary importance: as simulations predict a simple proportionality with the kick velocity (σ obs ≈ V k /3.3, Merritt et al. 2009), it is clear that a population of HCSCs observed in the halo of a galaxy would open the door to a direct determination of the kick velocity distribution, therefore constraining the merger history of the host galaxy, the models of galaxy-assembly, the simulations of merging BHs and the assumptions upon which they rest. From a determination of M and V k one could also derive γ, gaining insights on the distribution of stars in the nucleus at the time of the merger.
The predicted number of HCSCs bound to the Milky Way (MW) ranges from a few tens to a few thousands: O'Leary & Loeb (2009) estimated that a MW-like galaxy which undergoes a hierarchical assembly, with no major mergers since redshift z = 1, should retain in its halo hundreds of HCSCs bound to BHs with masses in the range 10 3 ≤ M • ≤ 10 5 M , and tens of clusters bound to BHs with M • 10 5 M . These HCSCs were ejected from the shallow potential-well of the building blocks of the main galaxy, and they were trapped in the region which collapsed to make the MW-like galaxy. Later, Rashkov & Madau (2014) used the cosmological simulation Via Lactea II (Diemand et al. 2008) to predict the properties of a population of relic intermediate-mass BHs (IMBHs) in the halo of a MW-type galaxy. These too are leftovers of the galaxy hierarchical assembly. They identified a population of "naked" IMBHs (their sub-haloes were destroyed during infall) and "clothed" IMBHs (residing in the nuclei of stripped galaxies). The naked BHs make up 40 to 50 per cent of the total population and are mostly located within 50 kpc of the halo center, where tidal stripping is more effective. These BHs are also associated with compact stellar clusters, not necessarily bound to recoiling BHs, but simple residuals of a stripped nuclear star-cluster. The total number of the relic population ranges between 70 and 2000, depending on the BH seeding scenario and the steepness of the M • -σ relation (Ferrarese & Merritt 2000;Gebhardt et al. 2000) adopted to populate the nuclei. They would be spatially-resolved and with apparent magnitude as bright as m V ≈ 16, in the scenario producing the most massive IMBHs.
Considering the challenges that an all-sky survey would imply to search for HCSCs in our own galaxy, Merritt et al. (2009) produced a prediction for nearby galaxy clusters, where the limited angular extension would ease the task: they estimated that Virgo should contain one HCSC with apparent magnitude K ≤ 20 and up to 150 with K ≤ 26.
Observational evidence for the existence of HCSCs and recoiled SMBHs is scant. O'Leary & Loeb (2012) mined SDSS DR7 finding 100 HCSC candidates with photometric properties consistent with theoretical models, however, to date none of these has been confirmed as a genuine HCSC. Postman et al. (2012) argued that the SMBH of the Abell 2261 brightest cluster galaxy was ejected via gravitational recoil, and they identified four HCSC candidates located in the vicinity of the galactic nucleus -possible carriers of the putative ejected SMBH; later, Burke-Spolaor et al. (2017) showed that two of them are consistent with either dwarf galaxies or stripped nuclei, while the nature of the other two remains poorly constrained. The puzzling object HVGC-1 was tentatively interpreted as a hypervelocity globular cluster ejected from the Virgo cluster in the aftermath of a three-body interaction; the HCSC scenario was dismissed because of the low metallicity and the low velocity dispersion (σ ≤ 80 km s −1 , too low when compared with the putative kick velocity, V k ≥ 1025 km s −1 , Caldwell et al. 2014). Boubert et al. (2019) considered the possibility that the fastest star in Gaia DR2 is bound to a recoiling IMBH, but they favoured an interpretation where the measurement is spurious. Several other sources have been proposed as candidate recoiling SMBHs (presumably clothed with an HCSC), but alternative interpretations remain often equally viable and difficult to rule out conclusively (e.g. Bonning et al. 2007;Shields et al. 2009;Batcheldor et al. 2010;Robinson et al. 2010;Jonker et al. 2010;Civano et al. 2010;Tsalmantza et al. 2011;Eracleous et al. 2012;Koss et al. 2014;Lena et al. 2014;Menezes et al. 2014;Markakis et al. 2015;Chiaberge et al. 2017;Makarov et al. 2017;Kalfountzou et al. 2017;López-Navas & Prieto 2018).
With the aim of facilitating the task of their identification, we present spectroscopic and photometric renditions of HCSCs bound to recoiling BHs with a mass of 10 5 M . Photometric renditions are built upon the dynamical simulations of Merritt et al. (2009) andLoeb (2012). The effects of a population of blue stragglers are accounted for, in both photometry and spectroscopy. Methods and results are presented in Sec.2, where we provide details on spectroscopic simulations, on the derivation of colours for a number of publicly available datasets, and on the rendition of the cluster morphology, which includes the effects of kick velocity and dynamical ageing. In Sec.3 we discuss the derived properties of HCSCs, along with search strategies and challenges in their identification. We sum up and conclude in Sec.4.

METHODS AND RESULTS
In this section we provide details on the methods and tools used to simulate spectra, colours, and morphology of HC-SCs. Results are also shown.

Synthetic spectra
We simulated a set of hypothetical HCSC spectra for clusters with about 10000 stars -the number of stars expected for a cluster bound to a 10 5 M black hole ejected at low velocity (v k = 150 km s −1 ), and relatively young (time since the kick τ k < 10 7 yr). For each cluster we assumed a single stellar population and a single metallicity. However, the grid of models presented here allowed to explore the effects of different metallicities and ages of the stellar population on the cluster properties.
To simulate the integrated spectra of HCSCs we generated atmospheric models for the individual stars which make up the cluster, we derived the corresponding spectra via a spectral synthesis software, and we co-added the individual spectra to produce an integrated spectrum for the HCSC as a whole. Additional details on the steps of this process are given below.
As a starting point, we created a Hertzsprung-Russell diagram to represent every evolutionary stage present in the chosen stellar population. These diagrams were produced using the web interface cmd 1 v3.3 in conjunction with the "PAdova and TRieste Stellar Evolutionary Code" (parsec v1.2S, Bressan et al. 2012;Tang et al. 2014;Chen et al. 2015); the evolutionary tracks provided us with the physical parameters for each of the stars present in the star cluster.
In the cmd interface we chose the YBC version of the bolometric correction (Chen et al. 2019), and we included the effects of circumstellar dust adopting the following composition: 60 per cent silicate plus 40 per cent Aluminum Oxide for M stars, and 85 per cent Amorphous Carbon plus 15 per cent Silicon Carbide for C stars (Groenewegen 2006).
To generate the isochrones we adopted a Kroupa Initial Mass Function following a two-part power law (Kroupa 2001(Kroupa , 2002 and we generated a grid of models with metallicities Z = 0.0002, 0.002, 0.02 (solar), 0.03, 0.07, and ages of the stellar population τ = 1, 7, and 13 Gyr. To extract stellar parameters we inverted the cumulative mass function generated via cmd, extracting randomly the parameters until we reached the expected total mass bound to the recoiling BH. The stellar parameters were then used to create a series of atmospheric models using atlas9 (Kurucz 1970). atlas9 is a local thermodynamic equilibrium one-dimensional planeparallel atmospheric modeling software which uses opacity distribution functions to reduce the computational time. These atlas9 atmospheres are used to generate synthetic spectra for each stellar evolutionary stage using synthe, a suite of programs requiring input model atmospheres, chemical abundances, and a list of atomic and molecular species (Kurucz & Furenlid 1979;Kurucz & Avrett 1981); we adopted the atomic and molecular lines lists provided in the Castelli website 2 . Finally, individual stellar spectra were co-added to create a synthetic integrated-light spectrum for each of the star clusters of a given age and metallicity. The synthetic spectra were created with a wavelength coverage of 3000-24000Å at high resolutions (R∼ 500, 000), and then degraded to the desired velocity dispersions.
Blue stragglers were also included in the stellar population of the cluster with the following approach: for each isochrone of a given age, we compiled a separate set of theoretical models with ages in the range of 10-90 per cent the age of the isochrone in question. We then randomly extracted stars with masses 1 -2 times the mass of the main sequence turn-off, with the number of blue stragglers satisfying the relation by Xin et al. (2011) for Galactic open clusters: where N 2 is the number of stars within two magnitudes below the main sequence turn-off. With this approach the fraction of blue stragglers is approximately 0.3, 1, and 2 per cent in clusters with a stellar population of age τ = 1, 7, and 13 Gyr, respectively. Here we note that this approach naturally yields a population of yellow and red stragglers, i.e. evolved blue stragglers which left the main sequence (e.g. Kaluzny 2003;Leiner et al. 2016).
The resulting spectra are shown in Fig.D1 for a range of metallicities and ages, along with the best-matching spectra from the Pickles Atlas (Pickles 1998).
Transmission curves were taken from the Spanish Virtual Observatory Filter Profile Service 3 . Best-fitting curves to the predicted colour-colour loci are given in Appendix A.
skymaker is a simulator of astronomical images; it takes as input a parameter file and a catalog. The former specifies a number of parameters, such as seeing, pixel scale, telescope details, and central wavelength, among many others. The latter is a catalog of coordinates, magnitudes, and a code to discern between stars and galaxies. If a point-spread function (PSF) is not provided, then the software derives one taking into account blurring due to the atmosphere, to the telescope motion, diffraction and aberration features, effects due to optical diffusion (i.e. scattered light), and intra-pixel response (i.e. the sensitivity variation within a pixel). Once the objects in the catalog are rendered, a sky background is added along with Poissonian and Gaussian noise, representing photon-noise and noise due to electronics, respectively.
To create the input catalog for skymaker we adopted the stellar positions derived via N-body simulations by Merritt et al. (2009), Fig.D2. The simulation was realised with N = 1.5 × 10 4 equal-mass particles and assuming an initial power-law density profile ρ ∝ r −7/4 . An instantaneous kick of magnitude V k was delivered to the cluster along the −X direction at t = 0. We used the results obtained at t = 100 (physical units can be obtained by multiplying the N-body time by GM • /V 3 k , with G the gravitational constant, M • the black hole mass, and V k the kick velocity). Stellar locations were derived by inverting the cumulative mass profile derived from the N-body simulation. The process consisted of four steps: 1) we selected a random number uniformly distributed between zero and the total mass of the cluster; 2) by inverting the cumulative mass profile we selected a radial distance from the cluster centre; 3) we picked the N-body point located at that radial distance; 4) the star coordinate was extracted from a Gaussian distribution centred at the location of the N-body point, and characterised by a standard deviation given by 10 per cent the value of its 3D distance from the cluster centre.
The number of stars in the cluster, right after the where γ, here assumed to be 7/4, as in the N-body simulation, is the slope of the stellar density profile before the kick (ρ ∝ r −γ ), and r • is the radius containing an integrated mass in stars equal to twice M • . Eq.2 was derived by combining eq. 5, 7, and 8 in Merritt et al. (2009). The scaling relation adopted for r • is: which we derived by fitting the data presented by Merritt et al. in their Fig.12 for "power-law" galaxies, i.e. galaxies with a deprojected inner density profile steeper than R −0.5 (e.g. Lauer et al. 1995, but see also Graham 2013). Although "core" galaxies (those with a central density profile shallower than R −0.5 ), are more likely to result from a history of mergers (e.g. Ebisuzaki et al. 1991;Makino & Ebisuzaki 1996;Merritt 2006), and therefore more likely to be the progenitors of HCSCs, we prefer to use the scaling relation derived for "power-law" galaxies as it is better constrained at the low-mass end of the SMBH mass distribution, which is the focus of this work. The difference between the extrapolated values of r • is, however, small: for M • = 10 5 M , we get r •,core ≈ 0.8 r •,power−law ; furthermore, if one would extend Merritt et al.'s Fig.12 to lower masses, such difference would turn out to be much smaller than the scatter in the data points.
To derive magnitudes we inverted the evolved integrated initial mass functions produced via parsec along with the isochrones (Sec. 2.1). Stellar masses and the corresponding magnitudes were then extracted and randomly assigned to the N-body points. In this process we ignored the effects of mass segregation as O'Leary & Loeb (2012) found only moderate evidence that this process was taking place in their models, even though the dynamical simulations included stars spanning a factor 10 in mass. Similar results on mass segregation were obtained for simulations of globular clusters with a nuclear BH: the BH quenches mass segregation by scattering sinking particles out of the core (e.g. Baumgardt et al. 2004;Gill et al. 2008).
We assumed the telescope and seeing parameters indicated in Appendix C, a distance to the cluster of 10 and 30 kpc, and a BH mass of 10 5 M . Finally, we point out that an accurate representation of the cluster morphology should take into account the cluster expansion. Dynamical studies conclude that the radius enclosing a fixed number of stars evolves as r ∝ t 2/3 , however the evolution is expected to deviate from this law when the number of stars in the cluster becomes small (e.g. O'Leary & Loeb 2012, and references therein). Furthermore, Merritt et al. (2009) warn that the true evolution of the cluster is likely affected by a number of phenomena which are still poorly understood and not fully implemented in dynamical simulations. For these reasons, as Merritt et al., we ignored the cluster expansion.  Figure 1. Colour-colour plots for different instruments and different bands. Circles: expected colours for the HCSCs with metallicity as indicated in the colour bar, and with circle size proportional to the age of the stellar population (1, 7, or 13 Gyr); dashed black contours: randomly selected galaxies; solid grey contours: randomly selected stars. Contour labels indicate the fraction of data outside the contour. Details on the origin and selection of the plotted data are given in Appendix B (continues on next page).

Effect of different kick velocities
As eq.2 shows, the stellar mass bound to the recoiling BH right after the kick is a function of the kick velocity. The effects of different kick velocities on a HCSC bound to a 10 5 M black hole are shown in Fig. 2, which simulates a 40 second Pan-STARRS-DR1 exposure in the r-band, and a 116 second Euclid /NISP exposure in the J -band (these are the nominal durations of single exposures for the 3π Pan-STARRS survey, e.g. Schlafly et al. 2012, and for the wide field survey planned for Euclid, e.g. Carry 2018).

Cluster evolution after the kick
After receiving the gravitational-wave kick, the cluster is believed to evolve via resonant relaxation (Merritt et al. 2009;O'Leary & Loeb 2012). More specifically, a fraction of the stars is lost via tidal disruption events, and a fraction is ejected because of large-angle scattering (Henon 1969;Lin & Tremaine 1980). Merritt et al. (2009) andLoeb (2012) performed dynamical simulations to quantify the effects due to these processes, finding that the rate of stellar loss depends on the mass of the BH, with clusters around more massive BHs evolving more slowly than clusters around lower mass BHs.
To implement the time evolution we used the results of the Fokker-Planck simulation of O'Leary & Loeb (2012) for a cluster bound to a 10 5 M black hole. We reproduced their function using a smoothly-broken power law: where the amplitude A was matched to the number of stars bound to the BH at the time of the kick (τ k ≤ 10 6 yr), and the remaining parameters being t b = 6.79 × 10 6 yr, α 1 = −0.12, α 2 = 0.34, and ∆ = 1.76. The resulting images, derived from the xy panel of Fig.D2, are presented in Fig.3a and 3b, for Pan-STARRS and Euclid /NISP, respectively.

DISCUSSION
To narrow down the number of HCSC candidates from the mare magnum of galaxies and stars sampled in a survey one can start with a selection based on colour. For unresolved candidates, this might be the only parameter available, followed by an estimate of the absolute magnitude, kinematics, and by the stellar velocity-dispersion, if spectroscopy is feasible. Resolved candidates will offer the possibility of further cuts based on the morphology. We discuss these aspects in the following sections.  Fig.1 shows that the predicted optical and NIR colours of the clusters are not peculiar with respect to the general population of stars and galaxies. This is not surprising, as we have not implemented (and we are not aware of) any peculiar process which could be at play in these clusters, and which could leave a distinctive imprint on their stellar population.

Colours
In the optical, the predicted colours follow the locus of stars and galaxies. It is only in the JHK s colour diagram that HCSCs show an offset from the peaks of the distributions for stars and galaxies. It is for this reason that adding a NIR colour-cut to a set of candidates selected on the basis of their optical colours decreases the sample size by a factor ten. Still, the number of candidates obtained via a search based solely on optical and NIR colours is too large to be constraining. Muñoz et al. (2014) showed that a uiK s colour-colour diagram allows for a clean separation between stars, globular clusters, and galaxies. Such colour-colour diagram was used by Caldwell et al. (2014) to investigate the nature of the hypervelocity cluster HVGC-1, finding that it falls in the region defined by the globular clusters. We computed the expected uiK s colours for our sample of simulated HCSCs (Fig.1, central panel, bottom row); unfortunately, we found a large scatter which hampers the diagnostic power of this diagram in the search for HCSCs.
It is clear that a colour selection alone will not yield a suitable sample of candidates, unless it turns out that some peculiar process takes place into the exotic environment of such clusters, leaving a strong signature on its spectrophotometric properties, for example an enhanced rate of stellar mergers leading to an excess of blue stragglers. Preliminary simulations show that increasing the number of blue stragglers affects mostly the UV colours of sub-solar metallicity clusters; instead, adopting the GALEX (Bianchi & GALEX Team 1999) transmission curves, HCSCs with metallicity Z ≥ 0.02 (solar and above) have always colours in the range 6 ≤ (FUV -NUV) ≤ 8, independently of the presence of blue stragglers. Evolved stragglers falling on the red-giant branch affect optical and NIR colours introducing a higher variance in the result of the simulation (e.g. making a young cluster with solar-metallicity appear similar to a super-solar metallicity cluster with an old stellar population).
We draw the reader's attention to the following three points. First, sub-solar metallicity clusters might be unusual: HCSCs, being dislodged galactic nuclei, will likely show super-solar metallicity, and, being residuals of the galaxyassembly process, their stellar population might also be old (e.g. if a HCSC was ejected at redshift z=1, and no star-   Fig.3a for Euclid/NISP. formation took place in it, then its stellar population should be at least 7 Gyr old). While an old stellar population with super-solar metallicity might seem an odd combination, we note that several-times super-solar solar metallicity has been inferred for the broad-line regions of high-redshift quasars, e.g. Juarez et al. 2009. The second point to note is that, to date, the only instrument performing sky-surveys and sensitive in the near UV is Gaia (the G BP passband covers the wavelength range 330-680 nm, e.g. Weiler 2018; Maíz Apellániz & Weiler 2018), with no dedicated facilities planned for the foreseeable future. Unfortunately, also our predicted Gaia colours mostly follow the stellar locus. Finally, a pilot search through public databases shows that only a tiny fraction (below 1 per cent) of candidate HCSCs with optical and infrared coverage posses also good quality UV data.
We note that our models assume a single stellar popula-tion (including blue stragglers) and a single metallicity. On the other hand, HCSCs might share some properties with ancient nuclear star clusters (NSCs): if a NSC was present at the time of the recoil, the HCSC would consist of a subset of stars from the inner region of the NSC. The study of these compact objects is still in its infancy, and it is limited to the local Universe, however, it is becoming clear that they are characterised by complex stellar populations (e.g. Mastrobuono-Battisti et al. 2019;Neumayer et al. 2020). While our individual models for HCSCs might be simplistic, the overall set explores a range of ages and metallicities, including some extreme cases. Clusters with more complex star-formation histories fall within the colour-colour loci defined by the simple models presented here, with the scatter in the distribution being inversely proportional to the number of stars bound to the BH. By comparison, Merritt et al. (2012) proposed two models, the first with constant star-formation rate for the 5 Gyr preceding the GW-kick, the second with a single stellar population formed at the time of the kick; furthermore, they explored three metallicity histories: a constant solar metallicity, a constant subsolar metallicity (Z = 2 × 10 −4 ), and the estimated time-evolution for the Galactic Centre. The ugr colours derived here for SDSS can be directly compared with those derived by O'Leary & Loeb. The models produce consistent results, with the larger scatter recovered here in the red-side of the colour diagram being ascribed to bright evolved blue stragglers.

(2009) assumed a single stellar population and a metallicity following a Gaussian distribution; O'Leary & Loeb
To summarise, the combination of parameters producing models which deviate from the bulk of the population of stars and galaxies corresponds to super-solar metallicities and old stellar ages. In light of the observations of supersolar metallicities in the nuclei of high-redshift quasars, this combination, although unusual, might not be unrealistic, allowing to select ancient HCSCs.
In the following section we show that an HCSC with an old stellar population resembles K and M giant stars, and we highlight the spectro-photometric differences between the two classes of objects.

Spectra
In Fig.D1 we presented the model spectra derived as explained in Sec.2.1. Our derivation of the spectra did not take into account the internal dynamics of the cluster and the resulting shape of the absorption line profiles; for this we refer the reader to Merritt et al. (2009), where the authors dealt with this aspect in great detail. Rather, we will compare the model spectra with the observed spectra of single stars to identify the most probable interlopers, and to identify any difference -besides a high velocity dispersion (i.e. strongly broadened absorption lines) -that would allow to distinguish stars from clusters.
We compared the models with the stellar spectra of the Pickles Atlas (Pickles 1998), a library of 131 stellar spectra at solar abundance, including all normal spectral types and luminosity classes. We matched the binning of our models and the Pickles spectra, and we selected those resembling the most to the HCSC models (i.e. those producing a smaller residual when subtracted from the model). In Table 1 we indicate, for each metallicity and stellar population age assumed for the cluster, the best-matching spectrum (multiple spectra are indicated whenever the residuals differ by less than 10 per cent). In Fig.D1 we show the simulated spectra, the best-matching spectra from the Pickles Atlas, and the residuals.
From Table 1 we can see that HCSCs with a young stellar population (τ ≈ 1 Gyr) and sub-solar metallicity resemble G and F stars, either on the main sequence, or on the giant and sub-giant phase. At the other end of the spectrum, clusters with an old stellar population (τ ≈ 13 Gyr) and super-solar metallicity resemble giant K and M stars. Overall, most of the times, the simulated spectra resemble those of giant K stars. It is, therefore, likely that a colour-colour selection alone would be heavily contaminated by such objects.
How could we distinguish an unresolved HCSC from a giant star? Any distance information would help constraining the absolute magnitude of the candidates, applying a further selection cut on the sample. Moreover, the plots presented in Fig.D1 show often, but not systematically, a blue excess in the spectra of HCSCs. This is likely due to the contribution of the hottest main sequence stars and the blue stragglers in the cluster. We note, again, that the addition of a population of blue/yellow/red-stragglers acts as a confounding ingredient, adding stochasticity to the appearance of the cluster: simulations which do not include stragglers produce spectra which age as expected, i.e. with the peak shifting towards longer wavelengths and emission in the red portion of the spectrum becoming more and more prominent as the age of the stellar population increases. The addition of stragglers, instead, does not simply contribute with a population of bluestragglers, adding emission to the blue spectral component. Rather, this population also contributes with a few bright giant stars giving a major contribution to the integrated spectrum, ageing, instead of rejuvenating, the appearance of the cluster. We also note that the adopted number of blue stragglers (0.3, 1, and 2 per cent for τ = 1, 7, 13 Gyr) could be a lower limit: as explained in Sec.2.1, this fraction was derived for Galactic open clusters, where the stellar density is much lower than the one expected for HCSCs. The high density in the nuclei of HCSCs might enhance the production rate of blue stragglers.
To summarise, the simulated spectra of HCSCs often resemble those of K-type giant stars, however, the presence of a blue excess could be used to distinguish an unresolved HCSC from a single star. Fig.2, 3a, and 3b show rendered versions of the N-body model presented in Fig.D2. For this model the time since the kick is t = 100 × GM • /V 3 k ; although this quantity amounts to 2700 yr for M • = 10 5 M and V k = 250 km s −1 , this state of the cluster should be representative for a full relaxation time, i.e. approximately 10 6 yr. After that time, we assume that the cluster starts loosing stars as described in Sec.2.3.2. Fig.2 shows the effect of increasing kick velocities for a cluster with solar metallicity (Z=0.02), a stellar population of intermediate age (τ * = 7 Gyr), and located at distances of 10 and 30 kpc. The first two columns show the rendering of a Pan-STARRS r-band image. At the resolution of Pan-STARRS (median seeing θ r ≈ 1. 2 and pixel size p = 0. 25) the cluster corresponding to V k = 150 and 250 km s −1 is (barely) resolved to a distance of 30 kpc, showing evidence for a low-density envelope and a denser, compact core. For V k = 500 km s −1 the cluster is featureless, at visual inspection, and barely resolved at best. With the adopted saturation level and exposure time (reflecting the duration of a single exposure in the 3π Pan-STARRS survey, Schlafly et al. 2012), the core is saturated when V k = 150 km s −1 and for V k = 250 km s −1 at d = 10 kpc.

Morphology of a resolved HCSC
The last two columns of Fig.2 show the rendering of a Euclid /NISP J -band image. As NISP will be diffractionlimited, with a pixel size of p ≈ 0. 3, it is not surprising to find a much sharper image, with the cluster being resolved also for V k = 500 km s −1 at d = 10 kpc. The core is saturated in all six renderings. It must be stressed, however, that the rendering depends on the background and the instrument PSF, two parameters which, at the time of writing, are not well-known.
An observable quantity that one could derive for a HCSC and use it to mine catalogs, similarly to O' Leary & Loeb (2012), is the Petrosian radius, i.e. the radius where the local surface brightness equals the average surface brightness within that radius (Petrosian 1976). However, this parameter is not ideal for clusters which are resolved and which include bright stars in their halo. Such clusters will be treated as a collection of individual objects in a survey catalog. Moreover, for the cluster as a whole, the condition defining the Petrosian radius might be satisfied at multiple radii because the local surface brightness would not be monotonic. Therefore, this parameter can not be used for most of the scenarios explored here.

Kick signatures
From Fig.D2 it is clear that at t = 100 the cluster core is still flattened, with the major-axis aligned with the kick direction. However, this important morphological feature (it carries information on the kick direction) is lost in the rendered images presented in Fig.2, 3a, and 3b. The reason resides in the combination of stellar density profile and resolution effects: first, the cluster is characterised by a steep density profile; i.e. already at r 1 × GM • /V 2 k (in N-body units), or r 0.02 pc (for M • = 10 5 M and V k = 150 km s −1 ) the stellar density is 10 2 − 10 5 times lower than it is at the very nucleus (the exact value depends on the adopted density distribution before the kick, see Fig.3 in Merritt et al.; we used simulations produced from an initial density ρ ∝ r −7/4 ). Second, when magnitudes and stellar evolutionary stages are assigned to the N-body points, a few of them correspond to the evolutionary stages "red giant branch" (RGB), "core helium burning" (CHEB), or "early asymptotic giant branch" (EAGB). These are all bright stars (-5 < M J < 3, -1 < M r < 3), and the great majority of them are packed in the central density cusp. As a result, such bright stars are not resolved in the mock observations, and the convolution with the PSF erases the photometric asymmetry which, naively, would be expected from the spatial distribution alone.
To constraint the kick direction one could, therefore, rely on the asymmetry of the extended envelope, which evaporates as time goes by. As it is evident in Fig.D2, at t = 100 GM • /V 3 k (or after one relaxation time) stars are still distributed asymmetrically in the envelope: defining the envelope as the region outside the elliptical bulge with semimajor axes 13 and 10 (in N-body units), about 40 per cent of the envelope stars are located within the quadrant defined by position angles in the range −45 : 45 deg, where PA = 0 corresponds to the positive side of the X-axis, or the counter-kick direction. However, Fig.3a, and 3b, which are renderings of the xy panel of Fig.D2, show that already at t ≥ 10 7 yr such asymmetry looses its prominence. Projection effects, and higher recoil velocities (implying a lower number of bound stars), would only weaken the asymmetry and the conspicuity of the halo as a whole.

Light profile
What is the light profile that one would recover from an image of a resolved HCSC? Similarly to O'Leary & Loeb (2012), we derived, for the clusters rendered in Fig.3a and 3b, the logarithmic slope of the cumulative light profile: Γ i = d ln(I i )/d ln(r) ≈ ln(I i+1 /I i )/ln(r i+1 /r i ) and the average intensity within annuli of radius r i = 0.23, 0.68, 1.03, 1. 76, 3.00, 4.63, 7.43, 11.42, 18.20, 28.20 arcsec (the same adopted to produce the SDSS and Pan-STARRS catalogs). The width of the annuli is δ = 1.25r i −0.8r i . The derived slopes are given in Table 2; plots of the intensity profiles are shown in Fig.4.
The clusters presented in Fig.3a and 3b are nearby and ejected with low velocity. Therefore, they are well-resolved and individual stars in the halo can be discerned. However, Fig.2 shows that clusters ejected at higher velocities and located further away would be almost featureless and with a smooth profile. In Fig.4 we show the light profiles for two of those clusters, making a comparison with a Plummer (1911) profile, which is typical for globular clusters. The comparison shows that the profile of a HCSC observed with a resolution of about 1 (as for PanSTARRS and other ground-based surveys) is similar to a Plummer profile. From the light profile alone, the HCSC can, therefore, be misclassified as a globular cluster. On the contrary, observations at higher resolution (such as those expected for Euclid ) reveal the presence of a nuclear cusp, which deviates from the Plummer profile.
In the barely-resolved scenario, a first step towards a better understanding of the object nature might be the inspection of the residuals obtained after fitting and subtracting a PSF. We used galfit (Peng et al. 2002(Peng et al. , 2010 to apply this approach on the rendered Pan-STARRS images for a cluster ejected with V k = 500 km s −1 and located at d = 10 kpc (bottom left panel in Fig.2, the adopted PSF was a star rendered with skymaker, using the same parameters adopted to render the cluster). The result is shown in Fig.5: the residuals show clear deviations from a pure point source, suggesting a more complex nature of the object. The resolved structure is due to the presence of an extended envelope, which might be the host of one or more RGB stars. Although a small fraction of the stars in the cluster are in such evolutionary phase (i.e. < 1 per cent), because of their absolute magnitude, which can be as low as M r = −1 (M J = −5), they have a considerable weight in determining the appearance of the cluster, both in the well-resolved scenario (where the cluster appears as an unresolved nucleus embedded into an extended envelope), and in the barely-resolved case (where a single off-center star might be responsible for an elongated extension, resembling a slightly-resolved binary star or a compact galaxy and a foreground bright star).
The presence of the extended envelope depends from (at least) three factors: the age of the cluster, the kick velocity, and the metallicity. Ageing of a cluster is due to the combined action of dynamical ageing and ageing of the stellar population. Dynamical ageing is more important for lowmass clusters; e.g. O'Leary & Loeb 2012 found that clusters bound to BHs with mass M • ≥ 10 7 M lose a very small fraction of their mass over 10 10 yr, instead, the number of stars bound to a 10 4 M BH decreases dramatically with time. A HCSC bound to a 10 5 M BH is predicted to loose almost 90 per cent of its initial stars, however, also for V k = 250 km s −1 the HCSC should retain almost 3000 stars at the time of the kick, which translates into a cluster consisting of about 300 stars after 10 10 yr. The effect of the stellar population ageing is clear, instead, in Fig.3a and 3b, where, in general, the cluster envelope fades away while stars grow old; however, depending on the initial age of the stellar population, the cluster might brighten up (e.g. because of stars going through the giant phase). The kick velocity, V k , has perhaps the most important effect on the presenze of an envelope, as shown in Fig.2, with its size being dramatically reduced for V k > 250 km s −1 . Finally, metallicity plays a mild role, with the number of bright stars decreasing at higher metallicity (high-metallicity stars loose more mass via stellar winds, e.g. Trani et al. 2014).

Search strategies and challenges
Searching for a subset of objects in a large data set requires the capability to remove false positives, and to identify the desired targets with the minimum amount of follow up observations.
The most distinctive characteristic of HCSCs is their high velocity-dispersion in conjunction with their moderate absolute magnitude; this combination of parameters separates them well from the population of globular clusters, ultra-compact dwarf galaxies, and elliptical galaxies (Merritt et al. 2009). However, the spectroscopic information needed to populate such a diagram is not readily available for a large data set, and it is expensive to acquire. More realistically, a search for HCSCs would start with the mining of large public databases, with the aim of narrowing down the sample, and collect, in progressive steps and for a small subset of prime candidates, the necessary information to constrain their nature.
As shown in Sec. 2.2 and 3.1, in a colour-colour diagram HCSCs fall near or at the maximum of the distribution of stars and galaxies. Therefore, a selection based on optical colours alone will be heavily contaminated by false positives, and will not produce a sample of candidates small enough to allow the collection of additional data without investing substantial effort. It is clear, however, that a database search would greatly benefit from a multi-wavelength ap- Table 2. Logarithmic slope of the cumulative surface brightness associated with the profiles presented in Fig.4. Assuming M • = 10 5 M , V k = 150 km s −1 , d = 10 kpc, unless otherwise specified. The range of values was derived from projections of the clusters on the xy, xz and yz planes. proach: probing a larger portion of the candidate's spectral energy distribution would set more stringent constraints. For example, our experiments show that combining optical and NIR constraints allows to reduce the sample size by a factor ten with respect to the selection based on optical colours alone. Alternatively, one could opt for a search of old highmetallicity clusters which, in colour-colour diagrams, lie outside the region containing the bulk of the population of stars and galaxies. When the candidates are not resolved, then kinematic and parallax information, in conjunction with absolute magnitude estimates will help constraining the nature of the candidates and remove a number of false positives: e.g. in a search targeting HCSCs in the Galaxy halo one could remove sources showing no evidence for proper motion, these are likely extragalactic sources; moreover, the proper motion of HCSCs should be peculiar with respect to the proper motion of other objects in the vicinity. When searching for unresolved clusters, many interlopers will be K-type giant stars, followed by G and M giants. Given a low-resolution broadband spectrum, accurately flux-calibrated, then a blue excess in the spectrum of the candidates, at λ < 5000 , might help distinguishing between single stars and composite objects.
For candidates which are barely resolved, then binarystars (either visual or physical) will be an important population of interlopers. In that case, a shift in the photometric centroid from single-epoch observations at multiple wave-   Fig.3a and 3b (first and second column) and for two clusters from Fig.2 (third column) for Pan-STARRS (top row) and Euclid/NISP (bottom row). The plot shows the average intensity of light within annuli. The adopted metallicity is indicated on top of the panels. The time since the kick (τ k ) is indicated in the legend in years for the clusters from Fig.3a and 3b. Symbols with a solid black edge indicate the presence of saturated pixels in the nucleus. The vertical dashed line, in the top row, marks the typical FWHM/2 of the Pan-STARRS seeing in the r-band. The luminosity of the older cluster with metallicity Z = 0.002 increases because of a few stars on the early asymptotic giant branch. For comparison, a Plummer (1911) profile, representative of globular clusters, is plotted as a dashed grey curve in the third column. We stress that the inner flattening visible in the first and second column is due to the presence of saturated pixels. lengths, and a periodic astrometric shift would allow to remove false positives. However, the orbital period of resolved binaries can be very long, precluding the detection of periodic astrometric shifts. When the candidates are resolved, additional constraints can be placed on the light profile and, if bright offcenter stars are present, then Gaia might also detect the relative motion of the stars and allow to set dynamical constraints.

SUMMARY AND CONCLUSIONS
The possibility that hundreds of HCSCs populate galactic haloes is alluring. Their discovery would show that supermassive black holes do merge, and that gravitational waves can indeed remove BHs from their birthplaces (i.e. galactic nuclei). Their characterisation (i.e. measurement of velocitydispersion, stellar mass, and effective radius) would allow to determine the distribution of GW recoils, would cast light on the assembly of the host galaxy, and on the distribution of stars in the nuclear environment at the time of merger. While several candidate HCSCs have been reported, none has been securely identified. If predictions are correct, they have already been imaged in existing surveys. Identifying them is, however, non-trivial. In an effort to further our understanding of these objects and to ease the task of their identification, we expanded on the existing literature by producing broad-band spectra and photometric renditions of hypothetical HCSCs bound to recoiling BHs of mass 10 5 M . We used the renditions to derive light profiles, and the broad-band spectra to compute the colours that HCSCs would display in a number of recent databases.
Photometric renditions were based on the dynamical simulations presented in Merritt et al. (2009) andLoeb (2012), which we used to implement the spatial distribution of stars and the cluster evaporation. We produced images of the clusters as they would appear in the 3π Pan-STARRS survey and in the Wide Field survey of the instrument NISP on board the forthcoming Euclid space telescope. Images and light-profiles were derived for a range of kick velocities, distances, metallicities, dynamical ages and ages of the stellar population, with the inclusion of blue, yellow, and red-stragglers. While clusters ejected at moderate velocities (250 km s −1 ) and located at a distance of a few tens of kpc are barely resolved by current surveys, they can be resolved by Euclid, with both the instruments NISP and VIS (the Visible Imager, an instrument with a better resolution than NISP). When observed at a resolution of about 1 , typical of ground-based surveys, HCSCs located a few tens of kpc away and ejected with velocities above 250 km s −1 would appear featureless and with light profiles resembling those of globular clusters. The inner cusp of HCSCs can be revealed, instead, with sub-arcsecond resolutions, such as those achievable with Euclid. The capability to resolve the light profile is important when implementing a search, as it allows to place additional constraints on the candidates (colours, we showed, have limited constraining power). It is, therefore, desirable to develop a library of models exploring the full parameter space of BH masses, stellar masses, and cluster age.
We used stellar evolutionary models to generate a set of broad-band synthetic spectra and the corresponding colours that HCSCs would display. While optical colours were previously derived by Merritt et al. (2009) andLoeb (2012), we computed colours for additional bands, and for specific surveys, including recent ones, such as Gaia and NGVS, among others. We covered, therefore, a larger portion of the spectrum, providing more stringent constraints for the selection of a sample. We show that most of the times, in a colour-colour diagram, HCSCs fall very close to the peak of the distribution of stars and galaxies. The exception are high-metallicity clusters (Z = 0.07) with an old stellar population (τ = 13 Gyr). While this combination is somewhat unusual, it might provide a good representation for old HCSCs: being dislodged galactic nuclei, leftovers of the galaxy-assembly process, it is not unreasonable to argue that they might resemble the nuclei of those high-metallicity quasars observed at high redshift. The fact that most HCSCs fall very near the peak of the distribution of stars and galaxies, in colour-colour diagrams, implies that a search based on colours alone will be heavily contaminated by false positives. The set of simulated spectra allowed us to make a direct comparison with a library of observed stellar spectra to identify the most likely interlopers, and we found that HCSCs resemble, most of the times, K-type giant stars. However, often the spectra of HCSCs show a blue excess with respect to those of single stars. A possible distinctive signature which requires, however, the availability of spectra with carefully calibrated fluxes.
The ever increasing availability of databases opens new opportunities to searching for HCSCs, and Euclid will soon perform an unprecedented optical and NIR survey covering 40 per cent of the extragalactic sky. Results of the present paper can be used to select candidates across multiple databases, and to gain insights on their nature. While we focussed on the properties of HCSCs bound to 10 5 M BHs (among the most massive and bright expected within the MW halo), it is desirable to produce a comprehensive library of models, exploring the full parameter space of BH masses (especially below 10 5 M ), stellar masses, and properties of the stellar population.

NOTES
A package for the reproduction of predicted colours and renditions is available on Zenodo (https://doi.org/10.5281/ zenodo.3763444).  Figure D1. HCSC model spectra with metallicity and stellar population age as indicated in the plot title, and best-matching stellar spectra from the Pickles Atlas. Residual (model-star) and fractional residual (model-star)/star are plotted in the bottom panels. Stellar names are indicated following the scheme xxy, with xx the spectral type, and y the luminosity class. The nomenclature rxxy and wxxy refers to metal-rich and metal-weak stars, respectively.