Formation of LISA Black Hole Binaries in Merging Dwarf Galaxies: The Imprint of Dark Matter

The gravitational-wave (GW) signal from merging intermediate-mass black holes (IMBHs) in the mass range 10⁴–10⁶ M_⊙ is one of the best anticipated targets of the Lisa Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017). Given their relatively low mass, such black holes (BHs) have been linked to globular clusters (within which they can merge with stellar-mass BHs and produce LISA events; e.g., Fragione et al. 2018a, 2018b), massive star clusters (where they can be detected via tidal disruption events; e.g., Lin et al. 2018), and dwarf galaxies (DGs). DGs are extremely numerous and experience on average three major mergers in their lifetime (Fakhouri et al. 2010), therefore possibly providing several LISA signals. Recently, observational work has revealed a population of IMBHs through the discovery of low-luminosity active galactic nuclei (AGN) in DGs both in the local Universe (e.g., Reines et al. 2013; Baldassare et al. 2015) and at high redshift (z ≲ 2.4; Mezcua et al. 2018). The best-known example with multi-wavelength observations is RGG 118, hosting an IMBH with a mass of ∼5 × 10⁴ M_⊙ (Baldassare et al. 2015).

Providing theoretical predictions for the expected merger rates of IMBHs is thus highly relevant to LISA. Recently, a number of studies have attempted to address the evolution of IMBHs in merging DGs utilizing cosmological simulations (Tremmel et al. 2015; Habouzit et al. 2017; Bellovary et al. 2018). However, due to limited resolution, such contributions are unable to resolve not only the formation of the IMBH binary, which would require pc-scale resolution (Mayer et al. 2007; Pfister et al. 2017), but also the preceding pairing phase in the merger remnant (Chapon et al. 2013; Pfister et al. 2017).

The dark matter (DM) distribution in the central regions of DGs may significantly affect the orbital decay of IMBHs and the efficiency of IMBH binary formation in interacting DGs. This is because dynamical friction (DF) depends on the background mass distribution (Colpi 2014) which, in the case of DGs, can be dominated by DM even at small radii. In particular, DF in constant-density cores is known to lead to the stalling of sinking perturbers (Goerdt et al. 2006; Cole et al. 2012; Petts et al. 2016; see also Di Cintio et al. 2017). Therefore, whether or not the DM in DGs is described according to the Navarro–Frenk–White (NFW; Navarro et al. 1996) profile, expected in cold-DM (CDM) cosmologies, or follows shallower mass distributions (e.g., Weinberg et al. 2015; Read et al. 2016) becomes particularly relevant. Shallow DM density profiles can arise either due to galaxy formation processes (e.g., baryonic outflows triggered by supernova explosions; e.g., Governato et al. 2010), or as a consequence of modifications of the underlying DM model (e.g., self-interacting DM, Spergel & Steinhardt 2000; or fuzzy DM, Hui et al. 2017). Currently, it is at least clear that a wide diversity of DM distributions exists in DGs (Oman et al. 2017), and that in some cases models with nearly constant-density DM cores appear to better reproduce observed galaxy rotation curves (e.g., Brooks et al. 2017).

There has been no quantitative work aimed at elucidating the effect of the DM density profile on the orbital decay of IMBHs and the efficiency of IMBH binary formation in merging DGs. Here we investigate this novel aspect of the co-evolution of galaxies and their central BHs via a series of high-resolution controlled merger simulations of DGs embedded in DM halos with different density distributions.

2.1. Initial Conditions

We employed the methods described in Kuijken & Dubinski (1995), Widrow & Dubinski (2005), and Widrow et al. (2008) to generate self-consistent N-body models of axisymmetric DGs consisting of exponential stellar disks, central BHs, and extended DM halos whose density profiles followed the general form (e.g., Łokas 2002)

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{\rho }_{{\rm{s}}}}{{(r/{r}_{{\rm{s}}})}^{\gamma }{(1+r/{r}_{{\rm{s}}})}^{3-\gamma }},\end{eqnarray} \tag{ 1 }$$

where γ and r_s denote the asymptotic inner slope and the scale radius of the profile, respectively. The characteristic inner density, ρ_s, depends on γ, the halo formation epoch, and the present-day values of the cosmological parameters (we adopt the concordance Λ-CDM cosmogony and z = 0).

To investigate the degree to which the inner DM density distribution affects the formation of BH binaries during DG mergers, we varied γ in three otherwise identically initialized DGs: γ = 1, 0.6, and 0.2 (see also Kazantzidis et al. 2013). The value of γ = 1 corresponds to the NFW profile, whereas γ = 0.6 and 0.2 indicate a mild density cusp and a nearly constant-density core, respectively. These adopted shallow inner slopes of γ < 1 are well motivated as they resemble those of both observed (e.g., Oh et al. 2015) and simulated (e.g., Governato et al. 2010; Tollet et al. 2016) DGs.

Each DG comprised an exponential stellar disk and a central BH with masses M_d = 4 × 10⁸ M_⊙ and M_BH = 10⁵ M_⊙, respectively. These values are consistent with the distribution of stellar and BH masses in the DG sample of Mezcua et al. (2018), where they present the largest IMBH sample beyond z ∼ 0. Moreover, all DG models consisted of a DM halo with a virial mass M_vir = 2 × 10¹⁰ M_⊙ (corresponding to a virial radius r_vir = 70 kpc) and a concentration parameter c_vir ≡ r_vir/r_s = 20 (comparable to the median concentration value for a z = 0 cosmological halo at this mass scale; e.g., Macciò et al. 2007). The choice of M_vir, near the upper limit of the values suggested by empirical models of the stellar mass–halo mass relation (Read et al. 2017), was dictated by our desire to employ a relatively low DM-particle-to-BH mass ratio, to prevent numerical two-body heating.

The vertical scale-height and central radial velocity dispersion of the stellar disks were equal to z_d = 0.2 R_d and σ_R0 = 20 km s⁻¹, respectively. Adopting a typical value for the halo spin parameter λ = 0.04 (e.g., Macciò et al. 2007), we assigned a disk radial scale-length R_d = 1.1 kpc (Mo et al. 1998),⁴ comparable to the effective radius of RGG 118 (Baldassare et al. 2017). The above parameters have been chosen according to typical values used in DG simulations (Kazantzidis et al. 2011), and in forthcoming work a wider range of parameters will be explored.

The resolution of each component (DM halo, stellar disk, and BH) of the models is listed in Table 1. We use a typical time-step of 1 Myr, within a hierarchical leap-frog time-stepping scheme that allows time-steps as short as 5 yr. It is important to stress that, with 1.2 × 10⁸ particles, our high-resolution simulations have a mass resolution nearly two orders of magnitude higher than in recently published cosmological simulations. Each BH was represented by a single particle, placed at the center of the DG initially at rest (Tremaine et al. 1994). We evolved all DGs in isolation to assess the adequacy of our numerical choices, confirming that the models are not affected by two-body relaxation and artificial heating of the disk particles via interactions with the more massive DM particles. These test simulations also established that the residual motion of the BHs around the centers of our DGs does not affect the interpretation of our results, as it is much smaller than the relative separation of the BHs in the cases when a hard binary does not form (see below).

Table 1.  Particle Specifications

Resolution	High	Low
NDM, N⋆, NBH	5 × 107, 107, 1	106, 5 × 105, 1
mDM, m⋆, MBH	400, 40, 105	2 × 104, 800, 105
DM, ⋆, BH	4, 1, 1	326, 109, 109

Note. From top to bottom: particle number, mass (in M_⊙), and gravitational softening (in pc) of DM, stellar, and BH particles per galaxy.

Download table as:  ASCIITypeset image

Figure 1 presents the density profiles of our DG models after relaxing them in isolation for 0.5 Gyr, corresponding to 10 dynamical times at 1 kpc.

Zoom In Zoom Out Reset image size

2.2. Simulations

Figure 2 presents the stellar surface-density maps of the high-resolution mergers, soon after a remnant has formed, highlighting a qualitative difference between the cuspy (γ = 1) and shallow (γ < 1) DGs. The former, having a deeper potential well, produce less-extended shells and tidal streams following the merger, resulting in a more compact disky remnant. Tidal tails are also sharper in this case, suggesting that deep photometric observations of the outskirts of interacting DGs might provide information about the underlying DM potential.

Zoom In Zoom Out Reset image size

The differences among our various experiments are even more striking when comparing the total volume–density profiles of the remnants (Figure 3). In the high-resolution cases, the profiles are nearly flat in the γ < 1 cases for r ≲ 0.2 kpc, whereas the remnant of the cuspy galaxies has a profile even steeper than in the initial conditions in the same region (compare to Figure 1). These differences arise at radii ≲1 kpc. Hence, not surprisingly, they are not seen in the low-resolution runs, which appear to have an almost identical slope down to three times the stellar softening.

Zoom In Zoom Out Reset image size

The structure of the merger remnant has striking implications on the decay of the BHs. In Figure 4, we show the BH-separation evolution for all of the simulations. At the beginning, the orbital history of the BHs is, by construction, very similar among all runs. However, the different initial γ implies a different enclosed mass, which causes the orbits to eventually differ: this can be clearly seen already at the second pericenter and, most importantly, at the respective formation times of the remnants, highlighted by the vertical lines. Once the BHs are in the merger remnant, stalling occurs at 50–100 pc in the cored models, while the BHs continue to sink efficiently by DF in the cuspy remnant. Such a remarkable difference is the key result of this Letter.

Zoom In Zoom Out Reset image size

We recall that the stalling of extended perturbers, such as globular clusters and galaxy satellites, has been widely documented for cored DM halos (e.g., Petts et al. 2016). However, here it is significant that the stalling radius is much larger than the distance required for other processes to take over DF and promote decay, as we explain below. In order to show this, we highlight the different orbital-evolution stages of the two BHs by recalling that the BHs will merge on a global timescale τ_Merg given by the sum of three distinct timescales:

$$\begin{eqnarray}&&{\tau }_{\mathrm{Merg}}={\tau }_{\mathrm{Rem}}+{\tau }_{\mathrm{HB}}+{\tau }_{\mathrm{GW}}.\end{eqnarray} \tag{ 2 }$$

For a fixed-mass ratio and orbital configuration (type of orbit, initial separation, and first pericentric separation), the time for the merger to form a remnant, τ_Rem, depends on the mass in the central regions of the galaxies. It takes longer to form a remnant in the shallow-slope cases than in the NFW case. To compute τ_Rem, we searched for the first occurrence when we do not see two distinct cores in the surface stellar-density maps (Figure 2), obtaining, for the high-resolution simulations, τ_Rem = 6.56, 7.07, and 7.43 Gyr for the γ = 1, 0.6, and 0.2 case, respectively.

The second timescale, τ_HB, which is also connected to the central mass distribution and therefore to the density slope of the remnant, is the time that it takes to form a hard BH binary inside the remnant and, in the absence of gas, is driven by three-body encounters with passing stars. Lastly, τ_GW is the time needed by the hard IMBH binary to coalesce via GW emission. The hard-binary separation is given by Merritt (2013)

$$\begin{eqnarray}&&{d}_{\mathrm{HB}}=\displaystyle \frac{{{GM}}_{\mathrm{BH}}}{8{\sigma }^{2}},\end{eqnarray} \tag{ 3 }$$

with σ denoting the stellar velocity dispersion at 10 pc (computed at τ_Rem).

Using Equation (3), we obtain d_HB = 0.06, 0.05, and 0.03 pc for the high-resolution cases γ = 0.2, 0.6, and 1, respectively. As the stalling radius is much larger than d_HB in the cored remnants, it follows that there is no alternative process to promote orbital decay in this case. Therefore, defining τ_HB by identifying the first time the BH separation reaches d_HB, it follows that it is indefinitely longer than the Hubble time in the low-γ cases, and so is τ_Merg as a result.

We note that the two high-resolution γ < 1 cases behave similarly, as expected from the total density profiles of the remnants (Figure 3). Likewise, small differences are seen among all the low-resolution runs (the BH separation oscillates in all cases around 1 kpc, roughly three times the stellar softening) despite the different DM profiles. This is again not surprising, given the small differences in the total matter profiles. Because the high-resolution runs show that the stalling radius is close to ≲0.1 kpc, it follows that this is not resolved in the low-resolution runs, being comparable to the stellar softening length. Second, we verified that the lower number of particles in the low-resolution runs might start to play a role in the late phase of orbital evolution, whereas relaxation timescales are of order 100 t_Hubble in the high-resolution runs, suggesting that our results are robust and reflect the collisionless nature of the systems. Taken together, artificial relaxation and softening considerations cast doubts on studying BH dynamics in cosmological simulations (including most zoom-in runs), as their resolution is usually much coarser than in our high-resolution runs.

The threshold in resolved dynamics induced by gravitational softening is the reason why, even in the cuspy high-resolution remnant, BHs eventually stop sinking before they can reach the hard-binary stage. However, given sufficient resolution, sinking should continue until the hard-binary stage is reached, likely within t_Hubble, given the steep orbital-decay curve seen in Figure 4.

In the cuspy case, it is conceivable to compute τ_GW after the hard-binary stage has been reached. As we do not use a direct N-body code, we can rely on Sesana & Khan (2015), who provide an estimate of the coalescence timescale at the onset of the hard-binary regime. Owing to the low nuclear stellar densities in our remnants, for our NFW merger we obtain τ_GW ≫ t_Hubble, regardless of the orbital eccentricity. However, it is noteworthy to mention that the prototypical DG hosting an IMBH, RGG 118, possesses a compact bulge-like or nuclear star cluster (NSC) component at its center (Baldassare et al. 2017). Moreover, observations of nearby bulgeless DGs with the Hubble Space Telescope (HST) find an association of NSCs with AGN, suggesting that a central density peak is a feature of these systems (e.g., Kormendy & Ho 2013).

Thus, in the cuspy case, we can consider the effect of such central component using for simplicity a Hernquist (1990) model with M_Bulge = 0.01 M_d and r_Bulge = 100 pc (consistent with the constraints in Baldassare et al. 2017), and find τ_GW ≪ 1 Gyr. The addition of such component (a bulge-like structure or a NSC) appears thus to be crucial for the coalescence of the BHs in the case in which a hard binary can form. Recalling that IMBH–IMBH mergers should be detectable at least out to z ∼ 10, hence essentially through the entire cosmic history, BHs inside DGs with cupsy halos should give rise to a significant LISA event rate at IMBH mass scales.

Providing quantitative estimates for the number of LISA GW signals from the inspiral and coalescence of IMBHs from DG mergers is extremely difficult, as they depend, among other factors, on the BH occupation fraction in DGs, the rate of DG mergers, and the number of DGs.

The AGN occupation fraction in DGs is low, perhaps reflecting the impact of supernova feedback in low-mass galaxies, which removes the gas needed to fuel the central BH.⁵ If gas accretion is inefficient, BH mergers then become relatively more important as a BH-mass growth channel relative to larger galaxies, wherein accretion is dominant.

Also, DGs are believed to experience on average three major mergers between z ∼ 12 and 0 (Fakhouri et al. 2010), in contrast to larger numbers when assuming, e.g., Milky Way-sized galaxies (which experience on average five major mergers and several more minor mergers).

While the potentially low BH occupation fraction and the relative rarity of DG merger events might seem detrimental to a high LISA-detection rate, this is outweighed by the fact that DGs with M_vir ∼ 10¹⁰ M_⊙ are ∼100 times more numerous at z = 0 relative to Milky Way-sized halos, and even more so at high redshift (e.g., Reed et al. 2007).

However, even with a large integrated number of DG mergers, there remains the question of if (and how fast) a BH binary can form during a DG encounter, which is the focus of this Letter.

Our main finding is that the inner DM slope is a fundamental property for predicting the formation (or lack thereof) of an IMBH binary. More specifically, in our high-resolution simulations, only mergers between cuspy DGs (γ = 1) favor the formation of a hard BH binary, whereas, in the low-γ mergers, the BHs stall at a separation of ∼0.1 kpc (Figure 4). As a consequence, the rate of LISA signals from the coalescence of BHs from DG mergers is a function of the fraction of NFW-like DGs. We stress here that, even though we have only performed mergers between DGs with identical γ, only one NFW-like DG in the pair should suffice, because the inner DM slope of a merger remnant is always close to the steepest slope of the two (e.g., Fulton & Barnes 2001; Dehnen 2005). Furthermore, in DM models other than CDM, such as self-interacting or fuzzy DM, in which the formation of constant-density DM cores is a consequence of the physics of the DM candidate, we expect a dearth of LISA events at IMBH scales relative to CDM, hinting at the exciting possibility of using LISA event rates to probe the nature of DM itself.

Our results depend sensitively on resolution: in the low-resolution simulations, all of the BHs stall, independent of γ. This is an important point to keep in mind when interpreting results from cosmological simulations, where the resolution is necessarily limited.

It is obviously very difficult to obtain reliable measurements of the inner DM profile in observed DGs. Moreover, recent hydrodynamic simulations in the CDM cosmogony have shown that a wide variety of inner DM slopes arises owing to baryonic-feedback effects (e.g., Tollet et al. 2016), in which case LISA event rates at IMBH scales will probe the physics of galaxy formation at DG scales. Before then, the different stellar distribution arising in remnants with different DM profiles highlighted in Figure 2 suggest that wide-field deep photometry is also potentially a probe of the underlying halo structure of DGs. We will investigate quantitatively the latter method in a forthcoming work.

We note that we have included neither a stellar bulge nor a gaseous component in our models. A stellar bulge would in principle increase the central enclosed mass, thereby decreasing the DF timescale (e.g., Souza Lima et al. 2017; Tamburello et al. 2017). Moreover, the presence of a bulge would also accelerate coalescence, once (if) a hard binary forms. It is, however, currently very difficult to constrain the structural parameters of large samples of DGs in order to provide a quantifiable effect (Mezcua et al. 2018).

The addition of a gas disk is indeed critical, as it is not clear if the gaseous component would decrease (e.g., Mayer et al. 2007) or increase (e.g., Fiacconi et al. 2013; Tamburello et al. 2017) the DF timescale, an issue that in DGs is complicated further by the important effect of supernova feedback on the interstellar medium (e.g., Habouzit et al. 2017).

Properly modeling the gaseous component, together with BH accretion and feedback, would be also important in a complementary way: the low-γ DG-merger BHs, if they were able to accrete gas, could potentially shine for a long time as dual AGN with separations of ∼0.1 kpc. It would therefore be possible to detect them as low-luminosity dual AGN (e.g., Capelo et al. 2017), as their distance would be large enough to be resolved by, e.g., HST (at low redshift). It is actually not excluded that such AGN pairs could be already present, undetected, in the existing samples (e.g., Mezcua et al. 2018).

Given the long stalling time of the BHs in these remnants, it is also possible to have a third body merging with the remnant (Amaro-Seoane et al. 2010), with the fate of the triple system depending on the mass of the third body (a possible binary-hardening (ejection), if the third body is less (more) massive). Given the low average number of mergers experienced by DGs, it is however unlikely that triple DG systems would be common. We defer the investigation of all of these issues to future work, together with a detailed analysis of the IMBH–IMBH event rates in the LISA band, as such an analysis would be beyond the scope of this Letter.

We thank Thanos Anestopoulos, Monica Colpi, Massimo Dotti, Zachary Schutte, Nathan Secrest, and Alberto Sesana for fruitful discussions. T.T., P.R.C., S.K., and L.M. acknowledge networking support, and S.K. also acknowledges financial support from the COST Action CA16104, GWverse, on "Gravitational Waves, Black Holes and Fundamental Physics." P.R.C. acknowledges support by the Tomalla foundation. L.M.W. was supported by the Natural Sciences and Engineering Research Council of Canada through the Discovery Grant Program.