Unveiling the Population of Wandering Black Holes via Electromagnetic Signatures

Central to this work is the fact that SMBH dynamical evolution can be accurately tracked down to sub-kpc scales in this simulation suite. This is done via a sub-grid correction that accounts for unresolved dynamical friction from stars and dark matter (Tremmel et al. 2015). This frictional force is estimated by assuming a locally isotropic velocity distribution and integrating Chandrasekhar's equation (Chandrasekhar 1943) from the 90° deflection radius, r₉₀, out to the gravitational softening length, _g , ⁶ of the SMBH. This results in an acceleration applied to each SMBH particle in the simulation given by

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{\mathrm{DF}}=-4\pi {G}^{2}\rho (v\lt {v}_{\mathrm{BH}})\mathrm{ln}{\rm{\Lambda }}\displaystyle \frac{{{\boldsymbol{v}}}_{\mathrm{BH}}}{{v}_{\mathrm{BH}}^{2}},\end{eqnarray} \tag{ 1 }$$

where $\mathrm{ln}{\rm{\Lambda }}=\mathrm{ln}(\tfrac{{\epsilon }_{g}}{{r}_{90}})$, v_BH is the SMBH's velocity relative to the local center of mass velocity of the closest 64 star and dark-matter particles, ρ is the mass density, and G is the gravitational constant. This correction, combined with the high dark-matter and stellar mass resolution in the simulations, results in realistic sinking timescales for off-center SMBHs (Tremmel et al. 2015, 2018b).

We adopt the SMBH merger criteria derived by Bellovary et al. (2010), in which two SMBHs merge if they approach within 0.7 kpc of one another (or, twice the gravitational softening length of 350 pc) and are also mutually bound ($\tfrac{1}{2}{\rm{\Delta }}{\boldsymbol{v}}\lt {\rm{\Delta }}{\boldsymbol{a}}\cdot {\rm{\Delta }}{\boldsymbol{r}}$). The limit of 0.7 kpc is justified because at separations smaller than this distance we no longer trust the simulation to correctly model the dynamics of the system, as the gravitational force between the black holes becomes increasingly unresolved. When a merger occurs, the masses of the two SMBHs are added together and the resulting SMBH is given a position and velocity such that momentum is conserved (neglecting potential recoil due to gravitational wave emission).

SMBHs are seeded in gas that has reached high densities (3 m_p cm⁻³, 15 times the threshold for star formation in Romulus) while being free of metals and still at relatively high temperatures (9500–10,000 K). This method selects gas with high Jeans mass that is also collapsing quicker than the typical timescale on which such a dense gas particle will form a star (∼10⁶ yr). This is meant to roughly approximate the formation processes suggested by massive initial seeding models including direct collapse (Lodato & Natarajan 2007; Alexander & Natarajan 2014; Natarajan 2021) and results in SMBHs that are seeded preferentially at z > 5 within low-mass galaxies (Tremmel et al. 2017). Ricarte et al. (2021) show that the resulting z = 0 SMBH occupation fraction matches estimates based on current observations down to dwarf-galaxy scales.

The initial mass of a SMBH is set to 10⁶ M_⊙, with each SMBH consuming nearby gas particles to account for their mass. This mass is somewhat higher than what is typically assumed even for the most massive direct-collapse formation scenarios, but it is required that the mass be at least a few times larger than the background dark matter and star particles to avoid spurious scattering events (Tremmel et al. 2015). It is possible that such large masses affect the growth history of these SMBH, but Ricarte et al. (2019) showed that the SMBH growth depends mostly on the galaxy properties, rather than the SMBH mass. It may also affect their dynamical evolution, but Tremmel et al. (2018a) showed that even with their relatively large masses, the wandering SMBH population has dynamical friction sinking timescales much longer than a Hubble time in z = 0, MW-mass halos.

A SMBH's accretion rate is estimated according to a Bondi–Hoyle–Lyttleton prescription (Bondi 1952), modified to account for angular momentum support. Specifically, the accretion rate is given by

$$\begin{eqnarray}{\dot{M}}_{\bullet }=\alpha (n)\left\{\begin{array}{ll}\displaystyle \frac{\pi {\left({{GM}}_{\bullet }\right)}^{2}\rho }{{\left({v}_{\mathrm{bulk}}^{2}+{c}_{s}^{2}\right)}^{3/2}} & \mathrm{if}\ {v}_{\mathrm{bulk}}\gt {v}_{\theta }\\ \displaystyle \frac{\pi {\left({{GM}}_{\bullet }\right)}^{2}\rho {c}_{s}}{{\left({v}_{\theta }^{2}+{c}_{s}^{2}\right)}^{2}} & \mathrm{if}\ {v}_{\mathrm{bulk}}\lt {v}_{\theta },\end{array}\right.\end{eqnarray} \tag{ 2 }$$

where G is the gravitational constant, M_• is the SMBH mass, ρ is the ambient mass density, c_s is the ambient sound speed, v_θ is the local rotational velocity of surrounding gas, and v_bulk is the bulk velocity relative to the SMBH. The coefficient α(n) provides a boost to the accretion rate given by

$$\begin{eqnarray}\alpha (n)=\left\{\begin{array}{ll}{\left(\displaystyle \frac{n}{{n}_{\mathrm{th},* }}\right)}^{2} & \mathrm{if}\,n\geqslant {n}_{\mathrm{th},* }\\ 1 & \mathrm{if}\,n\lt {n}_{\mathrm{th},* }\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where n_th,* is the star formation number density threshold, below which the simulation no longer resolves the multi-phase interstellar medium, as in Booth & Schaye (2009). During the accretion process, thermal energy is injected isotropically into the surrounding gas particles, assuming a radiative efficiency of 10% and a feedback coupling efficiency of 2%. We adopt the same radiative efficiency when computing bolometric luminosities, such that ${L}_{\mathrm{bol}}=0.1{\dot{M}}_{\bullet }{c}^{2}$.

We follow the post-processing performed in Ricarte et al. (2021), using the Amiga halo finder to identify structures (Knollmann & Knebe 2009), pynbody for particle data analysis (Pontzen et al. 2013), and tangos for the construction of a database of galaxy and SMBH properties (Pontzen & Tremmel 2018). As in this previous study, we designate a SMBH as a "wanderer" if it is further than 0.7 kpc (twice the gravitational softening length) from the center of its host halo, which is identified using a shrinking spheres method (Power et al. 2003). This distance threshold is used because the center of the galaxy/halo becomes ill defined within 0.7 kpc at the gravitational force resolution of Romulus. Further, we want to avoid selecting wandering SMBHs that may be about to merge with another SMBH (see Section 2.1 above).

Accreting, wandering black holes, like those shown in Figure 1, can manifest as HLXs. These are off-nuclear X-ray sources with X-ray luminosities above 10⁴¹ erg s⁻¹ (Kaaret et al. 2001; Matsumoto et al. 2001; Gao et al. 2003). As the most extreme tail of the ultraluminous X-ray (ULX) population, these sources are unlikely to be produced by stellar mass objects and may instead represent accreting wandering SMBHs (King & Dehnen 2005). Masses in the intermediate black hole mass range (M_BH < 10⁶ M_⊙) are also supported by X-ray spectral fitting, which in turn imply relatively low blackbody temperatures (Miller et al. 2003; Davis et al. 2011).

Barrows et al. (2019) assembled a sample of 169 HLXs by cross-matching galaxies with known redshifts in the Sloan Digital Sky Survey (SDSS) with the Chandra Source Catalog (Evans et al. 2010). This sample includes sources with hard X-ray luminosities in excess of 10⁴¹ erg s⁻¹ in the 2–10 keV band up to z ∼ 0.9. Barrows et al. (2019) conservatively keep only those X-ray sources that are offset by ≥5 times the uncertainty of their positions relative to their host galaxy centroids, which varies from source to source. Most HLXs in this sample exhibit offsets of tens of kpc.

In Figure 2, we compare the distribution of offsets obtained from this work with Romulus25. In black, we plot the distribution of offsets from Barrows et al. (2019), with each source weighted by the factor (1 − f), where f is their estimated contamination fraction of each source. The red histogram is the distribution of offsets in the Romulus25 simulation obtained when emulating their selection criteria. For this comparison, we first select galaxies with rest frame r < 22.5 for all snapshots with redshifts z < 0.9, to mimic the pre-selection of galaxies in the SDSS. Then, we search for every SMBH within the virial radius of each selected halo, even if it is within a satellite halo. If its bolometric luminosity averaged over the past 30 Myr is at least 10⁴² erg s⁻¹, (that is, if its X-ray luminosity is at least 10⁴¹ erg s⁻¹ assuming a bolometric correction of 10%), we consider it detectable and save its offset from the center of the halo. The three-dimensional spatial distribution of HLXs is then projected onto the sky using an Abel transform, and each redshift snapshot n in the simulation between 0.05 < z < 0.9 is volume-weighted by redshift. ⁷ These distributions are normalized to integrate to unity.

Zoom In Zoom Out Reset image size

Comparing the red Romulus25 predicted distribution to the observed HLX distribution shown in black, we do indeed reproduce a population of luminous wandering SMBHs with separations of tens of kpc, as observed in Barrows et al. (2019). The most extreme offset SMBHs typically reside in the outskirts of the halos of massive galaxies, with stellar masses >10^11.4 M_⊙. However, unlike the observations, we find that this distribution should continue to increase with decreasing projected separation. This implies the existence of a much larger population of moderately offset sources missed observationally by current campaigns (Barrows et al. 2019). This is most likely a selection effect due to limited astrometric precision, as discussed further in Stemo et al. (2020).

Using the Amiga halo finder, we can assess the fraction of detectable HLXs in the simulation that are "true" wanderers bereft of any resolved subhalo, as opposed to those that only appear as wanderers because they reside in nearby satellite galaxies too faint to have been detected. In the blue histogram of Figure 2, we repeat the analysis done to produce the red histogram, but omit any SMBHs that actually reside in a subhalo detected by our halo finder. This distribution declines more rapidly with radius, indicating that the most offset HLXs likely reside in faint satellite galaxies. We plot the fraction of wanderers lacking a subhalo in the simulation as a function of projected halo-centric distance in the top panel. We find that this fraction stays above 50% until a projected separation of 15 kpc.

For HLXs with subsequently identified optical counterparts, Barrows et al. (2019) also considered the SMBH mass (M_•) to stellar mass (M_*) relation for their sources. Lacking a more direct measurement of M_•, they estimated M_• by assuming that the HLXs accrete with an Eddington ratio of 0.24, motivated by the average X-ray spectral index observed and a relationship between spectral index and Eddington ratio (Greene & Ho 2007). We perform the same selection on our Romulus25 HLXs, which have an identifiable subhalo (those that are part of the inventory in the red histogram, but not the blue one, in Figure 2). In Figure 3, we plot the SMBH mass to stellar mass relation for the subset of the Barrows et al. (2019) sample with identified stellar counterparts as solid blue circles, where error bars have been omitted for clarity. We plot the SMBH mass as a function of host stellar mass in the Romulus25 simulation in orange, computed in a different way in each of the panels. In the left panel, we plot each SMBH's accreted mass ⁸ in the simulation, which has been found to trace the M_• − M_* relation in Romulus even below the seed mass of 10⁶ M_⊙ and can be interpreted as a proxy for the true black hole mass in dwarf galaxies (Ricarte et al. 2019). (See Appendix A and Figure 7 for a comparison of accreted and raw SMBH masses among our wanderers.) Barrows et al. (2019) found that their sources were broadly consistent with AGN in low-mass galaxies, the relationship found by Reines & Volonteri (2015), shown as a green dashed line. Similarly, we find that the accreted masses of Romulus25 HLXs are also consistent with typical galaxies of the same mass. Here, the appropriate comparison is actually the relationship found by Schramm & Silverman (2013) shown in red, which is the relation used to calibrate the SMBH accretion and feedback parameters in Romulus. This reproduction is not entirely trivial, as environmental processes could have potentially modified this relationship for galaxies in these subhalos.

Zoom In Zoom Out Reset image size

In the right panel, we emulate the Barrows et al. (2019) proxy by "estimating" masses from SMBH luminosities assuming an Eddington ratio of 0.24. We find that this would have caused us to greatly underestimate Romulus25 masses, because Romulus25 SMBHs tend to have much lower Eddington ratios than this adopted value (Ricarte et al. 2019). That is, Romulus HLXs tend to be fainter than those in Barrows et al. (2019) at a given SMBH mass.

The overall good agreement between Romulus predictions and the observed HLXs lends us confidence in the robustness of the simulated properties of the off-center SMBH population. In subsequent sections, we explore the key electromagnetic properties of the wandering population in Romulus.

Dual AGN are ubiquitous in the Romulus universe, at least at low luminosities where we can build meaningful statistics. In Figure 4, we plot the probability of there being a second AGN in a galaxy with a (bolometric) luminosity at least L₂ > RL₁, given that there is already a first AGN with a luminosity of at least L₁ > 10⁴² erg s⁻¹. We plot the luminosity of these wanderers as a function of host stellar mass across several redshifts. Two brightness ratios are provided, R = 1/2 in orange and R = 1/10 in purple. Error bars are estimated via bootstrapping in each stellar mass bin. Although a bolometric luminosity threshold of 10⁴² erg s⁻¹ is quite low, we are unfortunately limited by the relatively small (25 Mpc)³ volume of Romulus25, which also does not produce many high Eddington ratio systems (Ricarte et al. 2019).

Zoom In Zoom Out Reset image size

In Romulus25, a galaxy hosting an SMBH with L_bol > 10⁴² erg s⁻¹ at z = 0.05 has a roughly 10% chance of hosting a second SMBH shining with at least one-tenth the first's luminosity. However, we predict that dual AGN should be much more common at higher redshifts, motivating deep surveys and sensitive instruments capable of testing this prediction. At present, searches for spatially resolved dual AGN in the X-ray are confined to z ≲ 0.05 and R ≳ 0.01 (Foord et al. 2019, 2021). Yet at z = 4, the majority of galaxies in Romulus with one SMBH above this luminosity threshold host another with at least 10% of its luminosity. We have tested increasing the luminosity threshold to L₁ = 10⁴³ erg s⁻¹, and did not find any noticeable difference aside from poorer statistics. SMBH accretion rates are averaged over 30 Myr in Figure 4, but it is possible for AGN to have shorter duty cycles than this. To interpret these results assuming that the AGN only shine a fraction D of the time, L₁ should be multiplied by a factor 1/D to conserve the average accretion rate, but the dual AGN probability must be correspondingly divided by the same factor.

If wandering SMBHs prove to be difficult to find and resolve individually, they may likely manifest as a "halo" of emission or excess counts in stacked images. In Figure 5, we compute the average emission profile of wandering SMBHs that one would obtain by stacking halos of similar mass. Once again, we average wandering SMBH accretion rates over a 30 Myr timescale, and further assume that 10% of the radiated energy is emitted in the 2–10 keV X-ray band (e.g., Hopkins et al. 2007). We perform an Abel transform to project average spherical profiles onto the plane of the sky, and apply the appropriate cosmological factors due to the evolving luminosity and angular diameter distances with redshift to derive this estimate. Each color represents a different bin in halo mass, as indicated by the legend.

Zoom In Zoom Out Reset image size

At z = 0.05, halos in the 10^11–13 M_⊙ mass range exhibit average wandering light profiles orders of magnitude above the cosmic X-ray background, plotted in gray. The solid horizontal line represents the cosmic X-ray background levels in the 2.0–10.0 keV band (Cappelluti et al. 2017). Dwarf galaxy halos (purple) have too few wanderers with too little luminosity, while the galaxy-cluster halo in RomulusC (orange) exhibits suppressed emission, as most galaxies are quenched and quiescent in the cluster environment (Tremmel et al. 2019). Unfortunately, although the RomulusC cluster hosts 1613 wandering SMBHs, they accrete less efficiently in its hotter halo (see also Ricarte et al. 2021, for further discussion of this phenomenon).

As thin dotted curves, we plot the average flux from central SMBHs, blurred assuming a Gaussian PSF with a FWHM of 05, appropriate for the center of the Chandra field of view. Additionally, as thin dashed curves, we plot the expected contribution from X-ray binaries (XRBs), for which details are provided in Section 3.1. At low enough redshift, these faint, extended halos of wanderer emission are much larger than the Chandra PSF and may extend to ∼10'' above both the X-ray background and host galaxies' XRBs. However, wanderers become too faint to detect above X-ray background levels by z = 1. This signal could potentially be detected by stacking X-ray observations centered on similar mass galaxies, masking out satellites, and comparing the resulting profiles to expectations from the field.

A TDE flare is caused by the destruction and accretion of a star torn apart by tidal forces when it approaches too close to a SMBH (Rees 1988). Several tens of TDEs have been observationally identified, and soon orders of magnitude more are expected from ongoing and future surveys such as the Zwicky Transient Factory and Vera Rubin Observatory (Strubbe & Quataert 2009; van Velzen et al. 2011; Bricman & Gomboc 2020). The number density of SMBHs as a function of mass is the basic ingredient needed to compute theoretical estimates of the volumetric TDE rate, and this has typically been derived from scaling relations between central SMBH masses and host galaxy properties (Stone & Metzger 2016; Stone et al. 2020). In this section, we demonstrate that the majority of SMBHs with masses ≲10⁷ M_⊙ are in fact wanderers, missed entirely by this kind of accounting. If these wandering SMBHs can retain stellar counterparts (such as their nuclear star clusters (NSCs)) and disrupt stars at comparable rates, then the majority of TDEs originating from SMBHs with M_• < 10⁷ M_⊙ may be centrally offset.

In Figure 6, we plot the fraction of SMBHs of a given mass in Romulus25 that are centrals, for three different redshifts. We find that this fraction drops rapidly at low masses, such that at z = 0.05, only one in twelve 10⁶ M_⊙ SMBHs are centrally located. The central fraction decreases with decreasing redshift, as more minor mergers build up the wandering population. If (i) these statistics are representative of the real universe, (ii) wandering SMBHs disrupt stars at a comparable rate to centrals, and (iii) offset TDEs are equally identifiable, then the majority of TDEs due to SMBHs with M_• < 10⁷ M_⊙ should be offset from their galactic centers.

Zoom In Zoom Out Reset image size

This may require a significant upward revision of the theoretically estimated TDE rates, which we estimate here. Stone & Metzger (2016) expressed the theoretical TDE rate as a function of SMBH mass as $N{\left({M}_{\bullet }/{10}^{8}\ {M}_{\odot }\right)}^{B}$, where N = 2.9 × 10⁻⁵ yr⁻¹ and B = −0.404. Consequently, we can estimate this correction factor via $C={\sum }_{\mathrm{all}}{M}_{\bullet }^{B}/{\sum }_{\mathrm{cen}}{M}_{\bullet }^{B}$, representing a sum over all SMBHs in Romulus25 in the numerator and a sum restricted to central SMBHs in the denominator. We find that the volumetric TDE rate when accounting for wandering SMBHs should be revised upward by a factor of C = 8.0 at z = 0.05, owing mostly to 10⁶ M_⊙ SMBHs, to which current surveys are not sensitive. This may even be an underestimate, because the black hole occupation fraction in dwarf galaxies is highly uncertain and may well be unity (e.g., Mezcua et al. 2018; Baldassare et al. 2020). On the other hand, this potential enhancement is larger than predicted by previous studies of ultra-compact dwarf galaxies, potentially due to differences in the SMBH occupation fraction (Voggel et al. 2019). At present, there exists one compelling offset TDE candidate, which is located 12.5 kpc from the center of a lenticular galaxy and exhibits a L ∝ t^−5/3 light curve characteristic of TDEs (Lin et al. 2018). The occurrence rate of offset TDEs can place joint constraints on SMBH dynamics, the SMBH occupation fraction, and the NSC occupation fraction.