Hunting for Wandering Massive Black Holes

First, we discuss our fiducial case of A2, where a massive BH moves at a constant velocity with a Mach number ${ \mathcal M }=2$ into hot plasma that has a density fluctuation with a characteristic wavelength λ = 5 R_BHL. In Figure 2, we show the 2D snapshots of the accretion flow at the plane of z = 0 (i.e., perpendicular to the net angular momentum vector) at three different elapsed times of t/t_BHL = 2, 5, and 35. In the early stage (t = 2 t_BHL), the supersonic gas flow is attracted by the gravitational force of the BH and forms a bow shock with a symmetric structure in front of the BH. As the density fluctuations reach within the BH influence radius (t = 5 t_BHL; middle panels), two streams from y > 0 and y < 0 collide at y = 0 and dissipate the linear momentum parallel to the y-axis. Because of the density asymmetry, however, nonzero angular momentum is left behind the colliding flows, and thus the denser flow from 0 < y/R_BHL < 0.7 accretes onto the BH, forming spiral arms and shocks. In the late stage after several dynamical timescales (represented at t = 35 t_BHL; bottom panels), the laminar flow with a spiral structure turns chaotic and turbulent. Since the gas is adiabatically compressed owing to the lack of radiative cooling, thermal pressure is not negligible. Therefore, the turbulent flow becomes subsonic (${{ \mathcal M }}_{\mathrm{tur}}\approx 0.5;$ third column) and the rotational velocity is sub-Keplerian (v_ϕ/v_Kep ≈ 0.3; fifth column). In this turbulent stage, the BH is fed not only through the disk but also by free-falling gas with substantially small angular momentum.

Zoom In Zoom Out Reset image size

Figure 3 shows the time evolution of gas accretion rate $\dot{M}/{\dot{M}}_{\mathrm{BHL}}$ through the sink cell at r = r_in (top panel) and mean specific angular momentum to the z-direction ${\bar{j}}_{z}/{j}_{\mathrm{Kep}}$ of the accreted mass (bottom panel). Blue solid curves correspond to our fiducial case. The accretion rate rises up to the BHL rate by t = 2 t_BHL and drops to $\sim 0.1\,{\dot{M}}_{\mathrm{BHL}}$ when density bumps enter within the BH influence radius and supply angular momentum of ≳0.8 j_Kep into the accreting matter. At t > 10 t_BHL, mass accretion approaches a quasi-steady state at a mean rate of $\langle \dot{M}\rangle \approx 0.15\,{\dot{M}}_{\mathrm{BHL}}$, though the angular momentum of the accreting matter has a large fluctuation with a mean value of $\langle {\bar{j}}_{z}\rangle /{j}_{\mathrm{Kep}}\approx 0.25$.

Zoom In Zoom Out Reset image size

Figure 3 also shows the dependence of the accretion flow and its angular momentum on Mach number (${ \mathcal M }=0.5$, 1.0, and 2.0) and wavelength of the density fluctuation (λ/R_BHL = 2.5 and 5.0), respectively. For all the cases, the overall behavior of the accretion flow is qualitatively similar to that in our fiducial case: the accretion rate initially increases to $\sim {\dot{M}}_{\mathrm{BHL}}$ and decreases to a quasi-steady value after the density bumps carry angular momentum within R_BHL. In Figure 4, we also calculate the frequency distribution of ${\mathrm{log}}_{10}(\dot{M}/{\dot{M}}_{\mathrm{BHL}})$ and ${\bar{j}}_{z}/{j}_{\mathrm{Kep}}$ during the quasi-steady state.

Zoom In Zoom Out Reset image size

With a higher Mach number, the average accretion rate in the quasi-steady state tends to be lower: $\langle \dot{M}\rangle /{\dot{M}}_{\mathrm{BHL}}\simeq 0.25$, 0.20, and 0.15 in the simulations of A0, A1, and A2, respectively. The angular momentum of accreting matter weakly depends on the Mach number, and the peak value is kept at $\langle {\bar{j}}_{z}\rangle \approx 0.3\,{j}_{\mathrm{Kep}}$. Besides, as shown in Figure 4 (solid curves), the width of the distributions becomes wider as the Mach number increases. This indicates that the accretion flow becomes more unstable and turbulent for higher values of ${ \mathcal M }$. Note that since ${\dot{M}}_{\mathrm{BHL}}\propto {\left(1+{{ \mathcal M }}^{2}\right)}^{-3/2}$, the accretion rate is reduced by a factor of ≃13.3 from the A0 run (${ \mathcal M }=0.5$) to the A2 run (${ \mathcal M }=2$).

With a shorter wavelength of density fluctuation, the flow pattern becomes more complex, although the absolute values of accretion rates and angular momentum do not change significantly (dashed curves in Figures 3 and 4). In Figure 5, we show the distribution of the gas density and velocity vector at an elapsed time of t = 20 t_BHL for the A1 (left) and B1 (right) runs, respectively. When the half wavelength is sufficiently larger than R_BHL as shown in the left panel, the incoming stream from 0 ≤ y ≤ λ/2 supplies mass and angular momentum with j_z > 0 (i.e., the counterclockwise direction) within R_BHL. On the other hand, in the right panel, the incoming stream from −λ/2 ≤ y ≤ −λ carries a larger amount of angular momentum and flips the direction of angular momentum (see also the bottom panel of Figure 3 at t = 20 t_BHL). As a result of the flow collisions around ∼R_BHL, the accretion flow turns highly turbulent, and thus the angular momentum distribution becomes wider.

Zoom In Zoom Out Reset image size

Next, we describe the properties of the accretion flow onto a moving BH, considering the time-averaged profiles of physical quantities. In the following, we show time-averaged values over 10 ≤ t/t_BHL ≤ 40.

Figure 6 shows the radial structure of the angle-integrated mass inflow (dashed) and outflow (dotted) rates for the A0, A1, and A2 simulations. These rates are defined as

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{in}}=-{\int }_{0}^{2\pi }{\int }_{0}^{\pi }\langle \rho \cdot \min ({v}_{r},0)\rangle {r}^{2}\sin \theta d\theta d\phi ,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{out}}={\int }_{0}^{2\pi }{\int }_{0}^{\pi }\langle \rho \cdot \max ({v}_{r},0)\rangle {r}^{2}\sin \theta d\theta d\phi ,\end{eqnarray} \tag{ 9 }$$

where $\langle \cdot \cdot \cdot \rangle$ means the time-averaged value. We also define the net accretion rate by $\dot{M}\equiv {\dot{M}}_{\mathrm{in}}-{\dot{M}}_{\mathrm{out}}$ (solid). Note that both the inflow and outflow rates are proportional to the area (∝r²) at larger radii where a uniform medium moves with a constant velocity without being affected by the gravitational force of the BH.⁵ Within the BH influence radius (r < R_BHL), the mass inflow rate starts to deviate from ${\dot{M}}_{\mathrm{in}}\propto {r}^{2}$ and approaches ${\dot{M}}_{\mathrm{in}}\propto {r}^{1/2}$, while the outflow rate decreases toward the center. As a result, the net accretion rate is nearly constant, and the accretion system is in a quasi-steady state. The radial dependence of the mass inflow rate ${\dot{M}}_{\mathrm{in}}\propto {r}^{1/2}$ is consistent with the result of simulations where mass accretion with a broad range of angular momentum occurs (Ressler et al. 2018; Xu & Stone 2019). We note that this accretion solution is different from those of self-similar RIAF solutions for a static BH (see also discussion below): ${\dot{M}}_{\mathrm{in}}\propto {r}^{0}$ (ADAF; Narayan & Yi 1995b) and ${\dot{M}}_{\mathrm{in}}\propto r$ (CDAF; Quataert & Gruzinov 2000; Inayoshi et al. 2018b).

Zoom In Zoom Out Reset image size

In Figure 7, we present the angle-averaged radial profiles of the density, rotational velocity, and temperature for the six models. For all the cases, the density and temperature begin to increase toward the center within the BH influence radius (r ≲ 2 R_BHL), where the accretion flow forms a sub-Keplerian rotating disk with a mean velocity $| {v}_{\phi }| \approx (0.2\mbox{--}0.4)\,{v}_{\mathrm{Kep}}$. Since the flow is not fully supported by the centrifugal force, the time-averaged inflow velocity is comparable to ∼(0.3–0.5)v_Kep. As the inflow rate in the quasi-steady state is approximated as ${\dot{M}}_{\mathrm{in}}\propto {r}^{1/2}$, the density follows $\rho \propto {r}^{-1}$ (see the top panel of Figure 7). Since radiative cooling is neglected in our simulations, the accretion flow is adiabatically compressed by the gravity of the BH and the temperature increases to the center following T ∝ r⁻¹, as expected from energy conservation. Note that this treatment is valid only when the BH is embedded in a low-density diffuse plasma so that the radiative cooling time is longer than the dynamical timescale at r ≃ R_BHL (and the orbital timescale for wandering BHs at the outskirts of galaxies; see Section 4). In Figure 8, we show the time-averaged angular profiles at r = 0.2 R_BHL for the same physical quantities shown in Figure 7. Although the density and rotational velocity increase around the equatorial plane, the accretion flow is no longer a geometrically thin disk structure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The power-law density profile (ρ ∼ r⁻¹) is qualitatively different from those of known RIAFs: ρ ∼ r^−3/2 for ADAF solutions (Narayan & Yi 1995b) and ρ ∼ r^−1/2 for CDAF solutions (Quataert & Gruzinov 2000; Inayoshi et al. 2018b). The overall properties of the accretion flow are similar to those discussed by Ressler et al. (2018) and Xu & Stone (2019), where the angular momentum of accretion flows is widely distributed.

Figure 7 also shows the dependence of the physical quantities on the Mach number and wavelength of density fluctuation. While the density and temperature hardly depend on the choice of λ, the density decreases and temperature increases with higher values of ${ \mathcal M }$. The density reduction simply reflects the dependence of ${\dot{M}}_{\mathrm{in}}$ on the Mach number due to the input of different angular momentum within the BH influence radius, as shown in Figures 3, 4, and 6. We note that the ${ \mathcal M }$ dependence of temperature is not true, but is caused by the radius being normalized by the BHL radius. In adiabatic gas, the temperature is given by the virial temperature independent of ${ \mathcal M }:$ $T\propto {GM}/r\propto (1+{{ \mathcal M }}^{2}){\left(r/{R}_{\mathrm{BHL}}\right)}^{-1}$. The amplitude of the rotational velocity is a fraction of the Keplerian velocity within r ≲ 2 R_BHL, though the rotation direction is more time dependent for shorter wavelengths, as shown in Figure 4.

We note that our simulations do not treat an explicit viscosity. As discussed in previous studies (Igumenshchev & Abramowicz 2000; Igumenshchev et al. 2000, 2003; Narayan et al. 2000; Quataert & Gruzinov 2000), the angular momentum of the accretion flow can be transported by turbulence excited by colliding flows. To analyze the effect, we calculate the r − ϕ component of the mass-weighted Reynolds stress,

$$\begin{eqnarray}&&{\tau }_{r\phi }=\displaystyle \frac{1}{4\pi }\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\pi }{\int }}\langle \rho {v}_{r}^{{\prime} }{v}_{\phi }^{{\prime} }\rangle \sin \theta {\rm{d}}\theta {\rm{d}}\phi ,\end{eqnarray} \tag{ 10 }$$

where ${v}_{i}^{{\prime} }\equiv {v}_{i}-\langle {v}_{i}\rangle$. In Figure 9, we show the radial profile of the Reynolds stress normalized by ρv_Kep² for the three cases. The Reynolds stress increases with the Mach number because the flow is more turbulent, and for ${ \mathcal M }\gtrsim 1$ it is approximated by τ_rϕ ≃ Ar⁻², where A is positive. This positive value of τ_rϕ indicates that the turbulent motions transport angular momentum outward. By analogy with the standard α-viscosity model (Shakura & Sunyaev 1973), we define the effective viscous parameter by

$$\begin{eqnarray}&&\hat{\alpha }(r)\equiv \displaystyle \frac{\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\pi }{\int }}\langle \rho {v}_{r}^{{\prime} }{v}_{\phi }^{{\prime} }\rangle \sin \theta {\rm{d}}\theta {\rm{d}}\phi }{\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\pi }{\int }}\langle \rho {c}_{s}^{2}{\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{Kep}}\rangle \sin \theta {\rm{d}}\theta {\rm{d}}\phi },\end{eqnarray} \tag{ 11 }$$

to quantify the strength of turbulent viscosity. In our simulations, we obtain $\hat{\alpha }(r)\simeq 0.2$ within r ≃ R_BHL. Therefore, turbulence transports angular momentum effectively even without MHD effects. Recently, Ressler et al. (2020) found that MHD and pure-HD simulations show similar properties of wind-fed accretion flows onto a BH in a nuclear region. In their situation, similarly to our simulations, mass accretion is allowed owing to a wide distribution of angular momentum provided stellar winds, even absent much angular momentum transport led by the MRI.

Zoom In Zoom Out Reset image size

Because of limitations in computing time, we do not extend our computational domain down to the BH event horizon scale (r ∼ r_Sch). Instead, we conduct two additional simulations with different locations of the innermost grid, at r_in/R_BHL = 0.16 and 0.32. Figure 10 shows the radial profiles of time-averaged and angle-integrated mass inflow rate (dashed), outflow rate (dotted), and net accretion rate $\dot{M}$ for each value of r_in. Within the BH influence radius, the inflow rate dominates the outflow rate, and the net rate becomes constant for all the cases. The normalization of the net accretion rate nicely scales with ${\dot{M}}_{\mathrm{in}}(r={r}_{\mathrm{in}})\propto {r}_{\mathrm{in}}^{1/2}$. In Appendix A, we describe the physical reason why the inflow rate depends on r^1/2 with an analytical model.

Zoom In Zoom Out Reset image size

In order to check whether radiative cooling matters, we compare the heating timescale to the cooling timescale. Since $\rho \,\simeq \,{\rho }_{\infty }{\left(r/{R}_{\mathrm{BHL}}\right)}^{-1}$ and $T\simeq {T}_{\infty }{\left(r/{R}_{{\rm{B}}}\right)}^{-1}$, the timescale for free–free emission at the rate of ${Q}_{\mathrm{br}}^{-}\propto {\rho }^{2}{T}^{1/2}\propto {r}^{-5/2}$ is estimated as ${t}_{\mathrm{cool}}\propto \rho T/{Q}_{\mathrm{br}}^{-}\propto {r}^{1/2}$. Since the main heating source in a RIAF is viscous dissipation, the heating timescale is given by ${t}_{\mathrm{vis}}\,\simeq \,{\left(\gamma (\gamma -1)\hat{\alpha }{\rm{\Omega }}\right)}^{-1}\propto {r}^{3/2}$, where $\hat{\alpha }\,\simeq \,0.2$ and Ω ≃ 0.3 Ω_Kep for our case. Thus, the ratio of the two timescales is estimated as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{t}_{\mathrm{vis}}}{{t}_{\mathrm{cool}}} & \simeq & \,0.14\,{\left(\displaystyle \frac{\hat{\alpha }}{0.2}\right)}^{-1}\,\left(\displaystyle \frac{r}{{R}_{\mathrm{BHL}}}\right)\left(\displaystyle \frac{{\dot{m}}_{{\rm{B}}}}{{10}^{-3}}\right)\\ & & \times {\left(\displaystyle \frac{{T}_{\infty }}{{10}^{7}{\rm{K}}}\right)}^{-\tfrac{1}{2}}{\left(1+{{ \mathcal M }}^{2}\right)}^{-2}.\end{array}\end{eqnarray} \tag{ 12 }$$

Since the heating timescale is shorter than the cooling timescale everywhere within r ≃ R_BHL, radiative cooling does not play an important role in the accretion flow as long as ${\dot{m}}_{{\rm{B}}}\lesssim 7\times {10}^{-3}$.

The r_in dependence of the net accretion rate affects the actual BH feeding rate and radiative output from the nuclear disk at r ≪ r_in. Numerical simulations of RIAFs find that the positive gradient of the inflow rate (i.e., $s\equiv d\mathrm{ln}\dot{M}/d\mathrm{ln}r\gt 0$) ceases and the net accretion rate becomes constant within a transition radius of r_tr ≈ (30–100)r_Sch (e.g., Abramowicz et al. 2002; Narayan et al. 2012; Yuan et al. 2012; Sadowski et al. 2015). Assuming s = 1/2, the reduction factor of the net accretion rate is estimated as ∼(r_in/100 r_Sch)^s ≃ 5 for a RIAF onto a moving BH with ${c}_{\infty }=500\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and ${ \mathcal M }=2$. In Appendix B, we discuss how radiation spectra are modified by this effect.

We note that the similarity between MHD and HD simulations seen at larger scales would not hold all the way down to the event horizon scales. In the inner region (r ≪ r_in), since the adiabatic index of gas changes from 5/3 to 4/3 because of relativistic effects and cooling processes (synchrotron and/or thermal conduction), magnetic field would be dynamically more important as seen in MHD simulations with general relativistic effects. However, the estimation of the transition scale is beyond our scope in this paper.

In the framework of hierarchical structure formation in the ΛCDM model, lower-mass galaxies form first, and they subsequently merge to build larger objects. In this paradigm, massive elliptical galaxies in the local universe are expected to experience a large number of galaxy mergers in a Hubble time. As a natural result of multiple dry mergers at low redshifts (gas-rich mergers at high redshifts), binary SMBHs form at the galactic core, and some of them merge into a single SMBH through multibody BH interactions that likely eject the smallest BHs from the core (Ryu et al. 2018; Zivancev et al. 2020). Therefore, some ejected BHs, depending on the kick velocity, are still bound within the galactic halo and orbit at velocities of ${v}_{\infty }\sim {\sigma }_{\star }$, where σ_⋆ is the stellar velocity dispersion. When the orbiting BHs are fed with the diffuse gas of the surrounding host, they emit nonthermal radiation, as discussed below.

As an example of a massive elliptical galaxy, we consider M87. To model the properties of gas surrounding a wandering BH, we adopt the Chandra observational data from Russell et al. (2015): the electron density ${n}_{e}\approx 0.11\,{\mathrm{cm}}^{-3}$ ($\rho \approx 1.84\times {10}^{-25}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$) and temperature $T\approx 1.9\times {10}^{7}\,{\rm{K}}$ (${c}_{{\rm{s}}}\approx 6.5\times {10}^{2}\,\mathrm{km}\,{{\rm{s}}}^{-1}$) for gas at a distance of $r\approx 2\mbox{--}3\,\mathrm{kpc}$ from the center. Since the mass of the central SMBH is as high as $\simeq 6\times {10}^{9}\,{M}_{\odot }$ (Gebhardt et al. 2011; Event Horizon Telescope Collaboration et al. 2019), the masses of the wandering BHs would be in the range ${M}_{\bullet }\approx {10}^{7-8}\,{M}_{\odot }$, which corresponds to BH mass ratios of q ≈ 10⁻³ to 10⁻². These mass ratios are reasonable for massive ellipticals that have frequently experienced minor dry mergers (see Figure 1 in Ryu et al. 2018). The orbital velocity of the moving BH is estimated as ${v}_{\infty }\sim {\sigma }_{\star }\simeq 321\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (the stellar velocity dispersion is taken from Babyk et al. 2018), corresponding to ${ \mathcal M }\,\simeq \,0.5$. Since this estimation is somewhat uncertain and the result is sensitive to the choice of ${ \mathcal M }$, as shown below, we treat the Mach number as a free parameter in the range of $0.5\leqslant { \mathcal M }\leqslant 2$. For reference, for a BH with M_• = 3 × 10⁷ ${M}_{\odot }$ moving at a velocity of ${ \mathcal M }=1$, the BH feeding rate is approximated as $\sim 0.2\,{\dot{M}}_{\mathrm{BHL}}=2.2\times {10}^{-7}\,{\dot{M}}_{\mathrm{Edd}}$ for the A1 run.

Figure 11 presents the radiation spectra with different BH masses of M_• = 10⁷, 3 × 10⁷, and 10⁸ ${M}_{\odot }$ and Mach numbers of ${ \mathcal M }=0.5$, 1.0, and 2.0. We also overlay the sensitivity curve of ALMA, assuming a distance of 16.68 Mpc for M87 (Blakeslee et al. 2009). For all the cases, the radiation spectra have peaks in the millimeter band at ν_p ≃ 100 GHz, where the ALMA sensitivity is the highest. The peak luminosity increases and exceeds the ALMA sensitivity with higher BH masses and lower Mach numbers.

Zoom In Zoom Out Reset image size

In Figure 12, we show the 100 GHz continuum luminosity as a function of BH mass M_• for different Mach numbers. The two horizontal lines correspond to the detection limits for ALMA (solid) and ngVLA (dashed), respectively. This shows that wandering BHs with M_• ≳ 2 × 10⁷ ${M}_{\odot }$ and ${ \mathcal M }\,\lesssim \,1$ could be detectable with ALMA. The detectable BH mass is reduced by a factor of 2–3 with ngVLA, whose sensitivity is one order of magnitude higher than that of ALMA. Note that if those BHs are wandering at larger distances of ∼5–10 kpc from the galactic center, where the plasma density is lower, their luminosities decrease and thus the detection threshold for the BH mass increases by a factor of ∼2–4.

Zoom In Zoom Out Reset image size

We apply this argument to other nearby massive elliptical galaxies, assuming the existence of wandering BHs at their galaxy outskirts. Taking the observational data from Russell et al. (2013, 2015) and Inayoshi et al. (2020), we estimate the properties of gas surrounding those BHs and quantify their predicted bolometric luminosities and 100 GHz flux densities. The errors of density and temperature are given by the maximum and minimum values at distances of r ≃ 2–3 kpc from the centers. We assume the mass of the wandering BH to be 1% of the central SMBH mass, and we choose ${ \mathcal M }=0.5$ (note that σ_⋆/c_s ≃ 0.5–0.8 for most cases in our sample). As shown in Table 2, the bolometric luminosities produced from wandering BHs are on the order of $\simeq {10}^{35}-{10}^{36}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and the flux densities at 100 GHz are 10⁻¹ to 10³ μJy. We note that the ALMA sensitivity at 100 GHz is ∼9 μJy for 1 hr on-source integration. Therefore, BHs, if any, wandering at the galactic outskirts could be detectable in M87 and NGC 4472. With the capability of ALMA, only a few nearby (≲20 Mpc) ellipticals are interesting targets for hunting wandering BHs.

Table 2.  Luminosity of Wandering BHs

Name	D (Mpc)	$\mathrm{log}({M}_{\bullet }/{M}_{\odot })$	${n}_{{\rm{e}}}\,({\mathrm{cm}}^{-3})$	$T\,(\mathrm{keV})$	$\mathrm{log}\dot{m}$	$\mathrm{log}({L}_{\mathrm{bol}}/\mathrm{erg}\,{{\rm{s}}}^{-1})$	$\mathrm{log}({F}_{{\nu }_{{\rm{p}}}}/\mu {Jy})$
M87	16.68	9.789 ± 0.027	0.114 ± 0.016	1.650 ± 0.050	−5.918 ± 0.104	38.268 ± 0.220	2.685 ± 0.115
NGC 507	70.80	9.210 ± 0.160	0.029 ± 0.009	0.965 ± 0.015	−6.742 ± 0.288	36.132 ± 0.695	−0.344 ± 0.618
NGC 1316	20.95	8.230 ± 0.080	0.033 ± 0.007	0.620 ± 0.010	−7.385 ± 0.181	33.923 ± 0.437	−1.551 ± 0.404
NGC 4374	18.51	8.970 ± 0.050	0.022 ± 0.005	0.595 ± 0.025	−6.778 ± 0.157	35.703 ± 0.486	0.422 ± 0.328
NGC 4472	16.72	9.400 ± 0.100	0.029 ± 0.010	0.785 ± 0.005	−6.418 ± 0.233	36.944 ± 0.543	1.620 ± 0.393
NGC 4552	15.30	8.920 ± 0.110	0.018 ± 0.004	0.455 ± 0.035	−6.738 ± 0.257	35.863 ± 0.592	0.621 ± 0.469
NGC 4636	14.70	8.490 ± 0.080	0.028 ± 0.011	0.485 ± 0.015	−7.029 ± 0.244	34.879 ± 0.545	−0.336 ± 0.443
NGC 5044	31.20	8.710 ± 0.170	0.050 ± 0.009	0.645 ± 0.015	−6.748 ± 0.254	35.644 ± 0.658	−0.294 ± 0.534
NGC 5813	32.20	8.810 ± 0.110	0.042 ± 0.009	0.585 ± 0.015	−6.661 ± 0.208	35.856 ± 0.563	−0.085 ± 0.396
NGC 5846	24.90	8.820 ± 0.110	0.042 ± 0.009	0.625 ± 0.015	−6.689 ± 0.210	35.859 ± 0.515	0.128 ± 0.396

Note. Column (1): galaxy name. Column (2): distance. Column (3): mass of the central BH. Columns (4) and (5): electron density and temperature at ∼2–3 kpc. Column (6): accretion rate normalized by the Eddington rate (${\dot{m}}_{}\equiv \dot{M}/{\dot{M}}_{\mathrm{Edd}}$) by the scaling relation from the simulation A1, assuming that the wandering BH mass is 1% of the central BH mass and ${ \mathcal M }=0.5.$ Column (7): bolometric luminosity. Column (8): radiation flux density at ν_p = 100 GHz. The distance and central BH mass are taken from Inayoshi et al. (2020, and references therein). The electron density and temperature are taken from Russell et al. (2015) for M87 and from Russell et al. (2013) for the others.

Download table as:  ASCIITypeset image

Finally, we generalize this argument for early-type, gas-poor galaxies of several morphological types and give an estimate of the millimeter luminosity from wandering BHs as a function of the stellar velocity dispersion σ_⋆. To characterize the gas density and temperature of the ambient environment of the wandering BH, we approximate the density profile with an isothermal β-model

$$\begin{eqnarray}&&\rho =\displaystyle \frac{{\rho }_{0}}{{\left(1+{r}^{2}/{r}_{{\rm{c}}}^{2}\right)}^{1.5}},\end{eqnarray} \tag{ 13 }$$

where r_c is the core radius, and the core density and gas temperature are estimated with Equations (22) and (23) in Zivancev et al. (2020) (scaling relations fitted with data from Babyk et al. 2018) as

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{\rho }_{0}}{{\rm{g}}\,{\mathrm{cm}}^{-3}}\right)=0.6\,\mathrm{log}\left(\displaystyle \frac{{\sigma }_{\star }}{100\ \,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)-23.8,\end{eqnarray} \tag{ 14a }$$

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{T}{\mathrm{keV}}\right)=2.50\,\mathrm{log}\left(\displaystyle \frac{{\sigma }_{\star }}{100\ \,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)-1.06.\end{eqnarray} \tag{ 14b }$$

We estimate the mass of the central SMBH using the M_•−σ_⋆ relation (Kormendy & Ho 2013) and set the mass of the wandering BH to 1% of the nuclear SMBH. As a reference, the orbital distance of the wandering BH from the galactic center is set to r = 3 r_c, and its velocity relative to the surrounding hot gas is set to ${ \mathcal M }=0.5$. We estimate the relation between the luminosity at ν_p = 100 GHz and the central velocity dispersion as

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{\nu }_{{\rm{p}}}{L}_{{\nu }_{{\rm{p}}}}}{\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)=7.6\,\mathrm{log}\left(\displaystyle \frac{{\sigma }_{\star }}{100\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)+32.1.\end{eqnarray} \tag{ 15 }$$

For distances comparable to that of M87, galaxies with ${\sigma }_{\star }\gtrsim 320\,\mathrm{km}\,{{\rm{s}}}^{-1}$ yield ${\nu }_{{\rm{p}}}{L}_{{\nu }_{{\rm{p}}}}\gtrsim {10}^{36}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, which can be detected by ALMA.

The existence of intermediate-mass BHs (IMBHs; see a recent review by Greene et al. 2020) with M_• ≳ 10⁴ ${M}_{\odot }$ in our Galaxy has been argued based on observations of high-velocity compact clouds (Oka et al. 2017; Tsuboi et al. 2017; Ravi et al. 2018) and theoretical/numerical studies (Volonteri & Perna 2005; Bellovary et al. 2010; Tremmel et al. 2018b). Tremmel et al. (2018a) predict that Milky Way−size halos would host ∼10 IMBHs within their virial radii, and that they would be wandering within their host galaxies for several gigayears.

We apply the same exercise as in Section 4.1 for wandering BHs with kiloparsec-scale orbits within the Milky Way. To model the properties of the hot gas surrounding the Milky Way halo, we adopt the results of the Suzaku X-ray observations (Nakashima et al. 2018), which estimate a plasma temperature of T ≃ 3 × 10⁶ K and an emission measure of EM ≃ (0.6−16.4) × 10⁻³ cm⁻⁶ pc. Based on these results, we adopt ${n}_{{\rm{e}}}=0.01\,{\mathrm{cm}}^{-3}$ as the gas density around wandering BHs.⁶ In Figure 13, we show the radiation spectra of wandering BHs with M_• ≈ 3 × 10⁴–3 × 10⁵ ${M}_{\odot }$ located at ∼10 kpc from Earth. The spectra in the millimeter band extend to lower frequencies, where (ng)VLA has the highest sensitivity. We could detect IMBHs down to M_• ≳ 10⁵ ${M}_{\odot }$ for ${ \mathcal M }\lesssim 1$.

Zoom In Zoom Out Reset image size

There is additional indirect evidence of the existence of hot gas in the Milky Way halo at distances larger than ∼50 kpc, based on observations of the Local Group dwarf galaxies with gas removed by ram pressure stripping (Grcevich & Putman 2009) and absorption lines of high-velocity clouds associated with the Magellanic Stream that is close to pressure equilibrium with a hot plasma (Fox et al. 2005). Those observations suggest a lower density for the hot gas halo ($n\approx {10}^{-4}\,{\mathrm{cm}}^{-3}$), which, if true, would imply that wandering BHs in the Milky Way halo are too dim to be detected.

Observations have identified IMBHs in low-mass dwarf galaxies with confidence down to M_• ≈ 10⁵–10⁶ ${M}_{\odot }$, and more tentatively for M_• ≲ 10⁵ ${M}_{\odot }$ (Greene et al. 2020; see also Mezcua 2017; Mezcua & Domínguez Sánchez 2020). Cosmological simulations studying the occupation fraction of IMBHs in dwarf galaxies find that a significant fraction of them are not centrally located but wander within a few kiloparsecs from the galaxy centers (Bellovary et al. 2019). This is expected for dwarf galaxies because of their shallow gravitational potential wells and the longer dynamical friction timescale for the wandering BHs. Multibody BH interactions and GW recoils due to BH mergers further contribute to the off-nuclear population of IMBHs (Lousto et al. 2012; Bonetti et al. 2019).

Recent radio observations of dwarf galaxies with the VLA by Reines et al. (2020) reported a sample of wandering IMBH candidates that are significantly offset from the optical centers of the host galaxies. Based on an empirical scaling relation between BH mass and total stellar mass, these authors argue that the candidate wandering BHs might have masses in the range M_• ≈ 10⁴–10⁶ ${M}_{\odot }$. With radio luminosities of ∼10²⁰–10²² W Hz⁻¹ at ν₀ = 9 GHz, the sources are radiating at $\gtrsim {\nu }_{0}{L}_{{\nu }_{0}}\simeq {10}^{-6}\,{L}_{\mathrm{Edd}}$. Based on the radiative efficiency model for RIAFs, this level of (bolometric) luminosity can be produced only when BHs accrete at relatively high accretion rates of $\dot{m}\approx {10}^{-4}$. However, at such a high accretion rate, radio synchrotron photons are heated to X-rays via inverse Compton scattering (Ryan et al. 2017).

The brightness of the radio emission could be explained by synchrotron radiation from nonthermal electrons accelerated in a relativistic jet instead of arising from a disk. Since the majority of the Reines et al. candidate wandering BHs are point-like sources at a resolution of ∼025 (which corresponds to a physical scale of ∼85 pc at the median distance of the sources), the jet age can be constrained to ≲10³ yr for an assumed jet propagation speed of ∼0.3c (e.g., Nagai et al. 2006; Orienti & Dallacasa 2008). The hypothesis of young jets seems consistent with their steep spectral indices (F_ν ∝ ν^−0.79), analogous to compact steep-spectrum sources (O'Dea 1998), although the spectral indices were estimated over a narrow frequency range (9–10.65 GHz).