SUPERMASSIVE BLACK HOLES IN A STAR-FORMING GASEOUS CIRCUMNUCLEAR DISK

To compute the evolution of the SMBHs embedded in the CND, we use the code Gadget-3 (Springel et al. 2001; Springel 2005). This code evolves the system by computing gravitational forces with a hierarchical tree algorithm, and it represents fluids by means of smoothed particle hydrodynamics (SPH; e.g., Monaghan 1992).

The SPH technique has difficulty resolving some hydrodynamical instabilities, such as Kelvin–Helmholtz or Rayleigh–Taylor instability, because it generates spurious forces on particles that are in regions with steep gradients of density (Agertz et al. 2007; McNally et al. 2012). A proper treatment of these instabilities is important for modeling processes such as star formation and the turbulence in the ISM. These instabilities are better resolved by codes that use a Eulerian grid-based technique, such as AMR; however, in these codes the orbits of massive particles experience strong and spurious perturbations, making massive particles follow unphysical orbits. For this reason, Eulerian grid codes are not suitable to study the orbital evolution of SMBHs embedded in a gaseous CND (Lupi et al. 2015 and the references therein) and we choose to use SPH rather than a Eulerian grid-based code.

To study how the formation of stars in a gaseous CND can affect the orbital evolution of a pair of SMBHs embedded in the CND, we created recipes for star formation, gas cooling, and heating due to supernovae, and implemented these recipes in the Gadget-3 code. Our recipes have certain key differences from those included in Gadget-3. In particular, in our recipes the gas can reach lower temperatures, star formation does not depend exclusively on gas density, and we do not assume the stellar feedback to be instantaneous. The recipes we implement resemble those of Stinson et al. (2006), which ere were used in the SPH code Gasoline (Wadsley et al. 2004), and those of Saitoh et al. (2008, 2010), which were used in simulations of isolated galaxies and galaxy mergers.

In our implementation, the cooling function is computed for an optically thin gas with solar metallicity. For temperatures between 10^4.8 and 10⁸ K we compute the cooling essentially as described by Katz et al. (1996), and for temperatures below 10^4.8 K we extend the cooling functions down to 10 K using the parameterization of (Gerritsen & Icke 1997).

To model star formation, we examine every SPH gas particle and select those that, for some time step ${\rm{\Delta }}t,$ satisfy the following three criteria:

$$\begin{eqnarray}&&{n}_{{\rm{H}}}\ \gt \ {n}_{\mathrm{min}}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\;\cdot \;{\boldsymbol{v}}\ \lt \ 0\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&T\ \lt \ {T}_{\mathrm{max}},\end{eqnarray} \tag{ 3 }$$

where is ${n}_{{\rm{H}}}$ the number density, T is the temperature, and ${\boldsymbol{v}}$ is the velocity. ${n}_{\mathrm{min}}$ and ${T}_{\mathrm{max}}$ are parameters fixed before each simulation (they are described in Section 2.2). The gas particles selected in this way are treated as candidates for star formation.

For each gas particle that satisfies these criteria, we compute the probability p of giving birth to a star of mass ${m}_{\star }$ (this mass is fixed throughout the simulation) as

$$\begin{eqnarray}&&p=\displaystyle \frac{{m}_{\mathrm{gas}}}{{m}_{\star }}\left(1-{e}^{-\displaystyle \frac{{C}_{\star }{\rm{\Delta }}t}{{t}_{\mathrm{form}}}}\right),\end{eqnarray} \tag{ 4 }$$

where ${m}_{\mathrm{gas}}$ is the mass of the gas particle, ${C}_{\star }$ is the constant that controls the star formation efficiency, and ${t}_{\mathrm{form}}$ is the star formation timescale, computed as the maximum of the cooling time ${t}_{\mathrm{cool}}=\rho \;\epsilon \;/{n}_{{\rm{H}}}\;{\rm{\Lambda }}(T)$ and the dynamical time ${t}_{\mathrm{dyn}}={({\rm{G}}\;\rho )}^{-1/2}.$ Finally, for each star formation-eligible gas particle, we draw a random number r between zero and one. If $r\lt p$ a new star particle of mass ${m}_{\star }$ is spawned and the new mass of the parent gas particle is computed as: ${m}_{\mathrm{gas},\mathrm{old}}-{m}_{\star }.$

We assume that every newly formed star particle is a Single Stellar Population (SSP) with a three-piece power-law initial mass function (IMF) as defined in Miller & Scalo (1979). From this IMF, and the parameterization of Raiteri et al. (1996) for stellar lifetimes, we compute the number of type II supernovae (SNe II) events that occur in each SSP as a function of time. Then we assume that for each of these SN II events the star particle deposits 10⁵¹ erg of energy into the closest 32 gas particles. We do not include any heating by SNe Ia, because the typical time of integration of our simulations is 30 Myr and the SN Ia timescale is on the order of 1 Gyr.

In order to model the pair of SMBHs embedded in a CND, we use the same initial conditions as Escala et al. (2005). In these initial conditions, the ratio between the mass of gas and the mass of stars in the CND is consistent with the mean value obtained from observations of the nuclear region of ultraluminous infrared galaxies (ULIRGs; Downes & Solomon 1998; Medling et al. 2014; Ueda et al. 2014). However, the radial extent of the CND is comparable to the radial extent of the smallest and more concentrated observed nuclear disks.

Initially, the CND follows a Mestel superficial density profile and has a mass ${M}_{\mathrm{disk}}=5\times {10}^{9}\;{M}_{\odot }.$ The radius of the CND is ${R}_{\mathrm{disc}}\;=\;400\;\mathrm{pc}$ and its thickness is ${H}_{\mathrm{disk}}\;=\;40\;\mathrm{pc}$. The CND is modeled with 235331 SPH particles, each with a gravitational softening of 4 pc and a mass of $2.13\times {10}^{4}\;{M}_{\odot }.$

The stellar component is initially distributed in a spherical bulge. It follows a Plummer density profile having a core radius of 200 pc and a mass within $r=\;400\;\mathrm{pc}$ ${M}_{\mathrm{bulge}}(r\lt 400\;\mathrm{pc})=5\;{M}_{\mathrm{disc}}.$ This stellar bulge is modeled with 10⁵ collisionless particles with a gravitational softening of 4 pc and a mass of $2.45\times {10}^{5}\;{M}_{\odot }.$

The black hole pair is modeled by two collisionless particles of mass ${M}_{\mathrm{BH}}=5\times {10}^{7}\;{M}_{\odot }.$ These particles are initially symmetric about the center of the disk, in circular orbits of radius 200 pc. The orbital plane is the plane of the disk.

For our recipe of star formation we need to set four free parameters: (1) the star formation efficiency ${C}_{\star }$, (2) the minimum density of gas required for star formation to occur ${n}_{\mathrm{min}}$, (3) the maximum temperature that a gas particle that can give birth to a star may have ${T}_{\mathrm{max}}$, and (4) the mass of the newly formed star ${m}_{\star }$.

In order to resolve the formation of a clumpy multiphase medium, we assume a high star formation density threshold ${n}_{\mathrm{min}}={10}^{4}\;{\mathrm{cm}}^{-3},$ and a low star formation temperature threshold ${T}_{\mathrm{max}}={10}^{3}$ K (Stinson et al. 2006; Ceverino & Klypin 2009). We set the mass of the star particles as half of the original mass of the gas particles. Only the parameter ${C}_{\star }$ is varied from between simulations to obtain different star formation rates. However, we restrict our selection of ${C}_{\star }$ to values that reproduce an average SFR that is consistent with the empirical Kennicutt–Schmidt relation (Kennicutt 1998; Escala 2015). In Figure 1, we show the surface gas density and SFR obtained for different values of ${C}_{\star },$ and compare this to the observed relation. The values of the parameter ${C}_{\star }$ that we use are 0.005, 0.015, 0.05, 0.15, and 0.5 (in Table 1 these are the runs C0005, C0015, C005, C015, and C05, respectively). In Figure 2 we also plot the time evolution of the gas mass of the disk and the mass on new stars for the five different values of ${C}_{\star }.$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  List of Simulations and Their Parameters

Label	Code	SF	${T}_{{\rm{f}}}$ [K]	${C}_{\star }$	${\epsilon }_{\mathrm{BH}}$ [pc]
C05	G3 Mod.	yes	25	0.5	4
C015	G3 Mod.	yes	25	0.15	4
C005	G3 Mod.	yes	25	0.05	4
C0015	G3 Mod.	yes	25	0.015	4
C0005	G3 Mod.	yes	25	0.005	4
C05_.04	G3 Mod.	yes	25	0.5	0.04
C015_.04	G3 Mod.	yes	25	0.15	0.04
C005_.04	G3 Mod.	yes	25	0.05	0.04
C0015_.04	G3 Mod.	yes	25	0.015	0.04
C0005_.04	G3 Mod.	yes	25	0.005	0.04
C05_.004	G3 Mod.	yes	25	0.5	0.004
C015_.004	G3 Mod.	yes	25	0.15	0.004
C005_.004	G3 Mod.	yes	25	0.05	0.004
C0015_.004	G3 Mod.	yes	25	0.015	0.004
C0005_.004	G3 Mod.	yes	25	0.005	0.004
SDH05	G3	yes	104	0.05	4
E05	G1	no	-	-	4

Download table as:  ASCIITypeset image

In order to see how the resolution of the gravitational forces affects the simulations, we ran ten simulations using lower gravitational softening for the black holes. In five of these simulations we set ${\epsilon }_{\mathrm{BH}}\;=\;0.04\;\mathrm{pc}$ and in the other five we set ${\epsilon }_{\mathrm{BH}}\;=\;0.004\;\mathrm{pc}$. We label these runs with the suffixes "_.04" and "_.004," respectively.

We also ran two additional simulations to compare orbital evolution using our star formation recipe to orbital evolution with different recipes. In the first additional simulation, we used the hybrid multiphase model for star formation implemented in Gadget-3 (Springel & Hernquist 2003; Springel et al. 2005a). This model assumes that the star formation is self-regulated and is parameterized only by the star formation timescale ${t}_{0}^{*}.$ We use the typical value for this parameter, ${t}_{0}^{*}\;=\;$ 2.1, which provides a good fit to the Kennicutt law (see Figure 1). The other simulation (run E05) uses more idealized gas physics. There is no star formation and the gas follows an adiabatic equation of state. This simulation corresponds to run A of Escala et al. (2005).

All of the simulations that we run are summarized in Table 1.

In Figure 3 we show the time evolution of the SMBHs' separation a for seven simulations. Five of these runs use our prescription for star formation, feedback, and cooling. The values of the star formation efficiency ${C}_{\star }$ (see Equation (4)) used in these runs are 0.005, 0.015, 0.05, 0.15, and 0.5, which in Figure 3 correspond to the purple, yellow, blue, green, and red continuous lines, respectively. The black continuous line corresponds to a simulation using Gadget-3's hybrid multiphase model for star formation (SDH05). The black dashed line corresponds to run A of Escala et al. (2005; E05), which uses idealized gas physics and does not consider star formation.

Zoom In Zoom Out Reset image size

From Figure 3 we see that the time it takes for the SMBHs to reach a separation comparable to the gravitational softening is in the range of 7–25 Myr. The fastest migration time corresponds to the SMBHs in the E05 run (black dashed line) and the slowest migration corresponds to the SMBHs in the SDH05 run (black continuous line). In the simulations with our prescriptions, the SMBHs reach this separation at a time between these two limits.

In Figure 4 we show a zoom of Figure 3 for times larger than 5 Myr. Here we see that the time it takes for the SMBHs to stabilize at a separation comparable with the gravitational softening is shorter for higher values of ${C}_{\star }.$ These times are 8, 7.9, 11.6, 16, and 17.2 Myr for runs of C0005, C0015, C005, C015, and C05, respectively (see table 1). Also we found that in the two simulations with the highest star formation efficiency, there are no large fluctuations in the SMBHs' separation after t = 8 Myr.

Zoom In Zoom Out Reset image size

From Figures 3 and 4, we concluded that the orbital decay of the SMBH pair occurs over a timescale that is at most ∼2.4 times longer in simulations using our recipes than in simulations using more idealized gas physics (E05).

It is important to note that in our simulations the orbital decay timescale varies by a factor of ∼2.1, while the star formation efficiency extends over two orders of magnitude. So the orbital decay timescale in our simulations has a weak dependence on star formation efficiency.

The effect produced by ${C}_{\star }$ on the orbital decay of the SMBH pair should be expected because the orbital evolution of the pair depends on the torques produced by the gas, and the structure of the gas is sculpted by the star formation, cooling, and SN II heating. In the next section we explore this issue, computing the gas distribution of the CND for different star formation efficiencies.

In Figure 5 we show the density distribution of the gas at t = 10 Myr for two different values of ${C}_{\star }.$ From this figure, we find that for higher values of ${C}_{\star }$ a greater proportion of the gas has a low density and the maximum gas density is lower. This is because with more vigorous star formation, high-density gas is more efficiently converted into stars. With a greater number of stars, there is a greater number of supernovae explosions, which heat the cold dense gas around the stars and drive the formation of a hot diffuse medium.

Zoom In Zoom Out Reset image size

The torques that the gas exerts on the SMBHs are density-dependent, and we expect that most of the torque experienced by the SMBHs is due to high-density gas. From Figure 5, we would expect that in simulations with smaller values of ${C}_{\star },$ where the portion of high-density gas is greater, the torque of the gas on the SMBHs is more intense and the in-spiral of the SMBHs is faster. However, Figure 6 shows that the torque produced by gas with a density greater than 10⁶ cm³ fluctuates repeatedly about zero, indicating that the high-density gas does not always extract angular momentum from the SMBHs, but also deposits angular momentum into the SMBHs. This is consistent with the large fluctuations that we observe on the separation of the SMBHs in Figure 4. For completeness, in Figure 7 we plot the torque produced by gas that is less dense than 10⁶ cm³. It is clear from this figure that the effect of the low-density gas on the orbits of the SMBHs is negligible.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The fluctuations of the torque exerted by the high-density gas on the SMBHs may be the result of the spatial distribution of this gas in the disk. To explore this, we show the distribution of gas with densities greater than 10⁶ cm⁻³ in the plane of the disk at t = 10 Myr for run C0005 in Figure 8 and for run C05 in Figure 9. From these figures, we find that gas with density greater than 10⁶ cm⁻³ is mainly concentrated in clumps, two of which surround the two black holes. In this clumpy CND, the SMBHs have encounters with these high-density gaseous clumps, and depending on the characteristics of the encounter, they gain or lose angular momentum. This is reflected in the fluctuations of the gravitational torque produced by the high-density gas on the SMBHs (Figure 6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From Figures 8 and 9, it seems that the main difference between run C0005 and run C05 is the number of clumps formed. This difference may be the cause of the different numbers of fluctuations in the SMBHs' separation seen in these runs, and the slightly different migration timescale of the SMBHs. Simulations with a lower SFR yield a greater number of gaseous clumps, so it is more probable for the SMBHs to closely interact with one or more clumps in these simulations.

In order to determine if the number of clumps depends on ${C}_{\star },$ we compute the number of clumps (${N}_{\mathrm{cl}}$) as the number of groups of gas particles that are gravitationally bound, with a central density greater than 10⁶ cm⁻³. We compute ${N}_{\mathrm{cl}}$ at three different times—5, 10, and 15 Myr—for each value of ${C}_{\star }$ explored. We also compute the mean density and the range of densities of these clumps. We summarize this information in Table 2.

Table 2.  Gas Clumps Number and Densities

${C}_{\star }$	${N}_{\mathrm{cl},{\rm{t}}1}$	${N}_{\mathrm{cl},{\rm{t}}2}$	${N}_{\mathrm{cl},{\rm{t}}3}$	${\rm{\Delta }}{n}_{\mathrm{cl}}\;({{\rm{cm}}}^{-3})$	$\bar{{n}_{\mathrm{cl}}}\;({{\rm{cm}}}^{-3})$
0.005	86	44	38	${10}^{6}-{10}^{8.7}$	${10}^{7.6}$
0.015	84	39	25	${10}^{6}-{10}^{9.3}$	${10}^{7.9}$
0.05	87	32	19	${10}^{6}-{10}^{8.6}$	${10}^{7.7}$
0.15	70	27	15	${10}^{6}-{10}^{8.8}$	${10}^{7.9}$
0.5	37	16	9	${10}^{6}-{10}^{8.0}$	${10}^{7.1}$

Note. ${N}_{\mathrm{cl},{\rm{t}}1}$ corresponds to the number of clumps at t1 = 5 Myr, ${N}_{\mathrm{cl},{\rm{t}}2}$ corresponds to the number of clumps at t2 = 10 Myr, and ${N}_{\mathrm{cl},{\rm{t}}3}$ corresponds to the number of clumps at t3 = 15 Myr. ${\rm{\Delta }}{n}_{\mathrm{cl}}$ corresponds to the range of densities of the clumps and $\bar{{n}_{\mathrm{cl}}}$ to the mean density of the clumps.

Download table as:  ASCIITypeset image

From Table 2 we found that the number of clumps decreases as ${C}_{\star }$ increases. This is consistent with Figure 5, where we found that higher SFRs correspond to a lower proportion of high-density gas. Therefore, in simulations with lower ${C}_{\star },$ the SMBHs are prone to interact with a higher number of gaseous clumps. This is reflected in the slight increase in fluctuations of both the separation and the gravitational torque experienced by the SMBHs that we observe in simulations with higher SFRs.

Having explained in this section why the orbital decay timescale changes with ${C}_{\star }$ in simulations using our recipe, in the next section we analyze why the orbital decay of the SMBHs in these simulations is faster than in simulation SDH05 and typically slower than in simulation E05.

As seen in Section 3.2, the stochastic fluctuations of the SMBHs' separation comes from gravitational interaction between the black holes and high-density gas clumps (${n}_{{\rm{h}}}\gt {10}^{6}\;{\mathrm{cm}}^{-3}$). However, these fluctuations are unlikely to occur in real CNDs because the intensity of the gravitational torques is density-dependent, and the densities of the gas clumps in our simulations are higher than the densities of gas clumps in observed CNDs. For example, if we compare our clumps' densities with the density of molecular clouds in isolated galaxies (Mathis 1990; Oka et al. 2001; Struve & Conway 2010) or ULIRGs (Downes & Solomon 1998; Schulz et al. 2007), which are typically $\sim {10}^{2}-{10}^{5}\;{\mathrm{cm}}^{-3}$, we find that the density of the clumps in our simulations is at least two orders of magnitude greater. Inside the molecular clouds, there are star-forming regions (cold or hot cores) with densities comparable to the mean density of the clumps in our simulations ($\langle {n}_{{\rm{h}}}\rangle \;\sim \;{10}^{7}\;{\mathrm{cm}}^{-3}$). However, they are one or two orders of magnitude less dense than roughly half of the gas clumps in our simulations (${n}_{{\rm{h}}}\;\sim \;{10}^{6}-{10}^{9.3}\;{\mathrm{cm}}^{-3}$). Even though one of the gas clumps that perturbs the SMBH orbits may have a density comparable with the density of these cores, the mass of these cores is on the order of $\sim {10}^{3}\;{M}_{\odot }$ (Garay et al. 2004; Muñoz et al. 2007), which is at least two orders of magnitude smaller than the mass of any of the clumps in our simulations.

The clumps in our simulations have these extremely high densities because our recipes do not account for all of the physical processes that sculpt the clumps. For example, in our simulations we use a cooling that assumes that the gas is optically thin while the clumps in our simulations are optically thick ($\tau \gg 1$). Also, there is a large amount of energy that we are not considering in our simulations, such as the energy of photons emitted by the stars and by the black hole's accretion disk. The effects of turbulence, which is typically damped in SPH codes, are also omitted. All these physical processes will prevent the further collapse of the clumps, promoting the formation of clumps with lower (and hence more realistic) densities.

The high densities of the gaseous clumps can affect the orbits of the SMBHs in a spurious way. The gravitational pull at the edge of these clumps can be greater than the gravitational pull at the "edge" of the SMBHs in our simulations, where we define the edge of an SMBH to be its gravitational softening. We define the effective mass density of our black holes (${\rho }_{\mathrm{BH}}$) to be the mass of the black hole (${M}_{\mathrm{BH}}$) enclosed within a sphere of radius equal to the gravitational softening of the black holes (${\epsilon }_{\mathrm{BH}}$), i.e., ${\rho }_{\mathrm{BH}}\;=\;3\;{M}_{\mathrm{BH}}/(4\pi {\epsilon }_{\mathrm{BH}}^{3}).$ We can compare ${\rho }_{\mathrm{BH}}$ with the mass density of the gaseous clumps (${\rho }_{\mathrm{cl}}\;=\;{m}_{{\rm{H}}}\;{n}_{\mathrm{cl}}$) to determine if, in our simulations, the clumps are more compact than the black holes. For a gravitational softening of 4 pc, we find that ${\rho }_{\mathrm{BH}}\;=\;$ 3.8 ${m}_{{\rm{H}}}\times {10}^{6}\;{\mathrm{cm}}^{-3}$. Hence, in our simulations the density of the gaseous clumps is greater than the effective mass density of the black holes. This implies that the gravitational pull at the edge of the black holes (${F}_{\mathrm{BH}}(\epsilon )\;\propto \;\epsilon {\rho }_{\mathrm{BH}}$) is typically smaller than the gravitational pull at the edge of the clumps (${F}_{\mathrm{cl}}({R}_{\mathrm{cl}})\;\propto \;{R}_{\mathrm{cl}}\;{\rho }_{\mathrm{cl}}\;=\;4-5\;\epsilon \;{\rho }_{\mathrm{cl}}$). Indeed, considering a gravitational softening of 4 pc, and the minimum and maximum densities of the gaseous clumps in our simulations, ${F}_{\mathrm{BH}}(\epsilon )/{F}_{\mathrm{cl}}({R}_{\mathrm{cl}})\;\sim$ 0.04–10. Therefore, a black hole sees the clumps as point masses, and it is not able to disrupt them when it has a close encounter with one of them, within a distance of ≲4 pc. Instead, in our simulations a close encounter within this distance will result in the scattering of the black hole.

If we decrease ${\epsilon }_{\mathrm{BH}},$ we would expect the higher resolution of the black hole's gravitational force to allow the black hole to disrupt a clump in a close encounter, where the minimum distance is ≲4 pc. To illustrate how decreasing ${\epsilon }_{\mathrm{BH}}$ affects the outcome of close encounters, we show how an encounter proceeds for different values of ${\epsilon }_{\mathrm{BH}}$ in the bottom two rows of Figure 11. Snapshots that are farther to the right in this figure correspond to later times. The first row of this figure shows a different encounter—one that illustrates scattering of the SMBH. In the middle row (run C0015_.004), we see that the SMBH's gravitational potential disrupts the clump, and the SMBH does not scatter significantly. On the other hand, in the lower row (run C0015) where ${\epsilon }_{\mathrm{BH}}$ is larger, the black hole is not able to disrupt the clump. The orbit of the SMBH is perturbed and the gaseous clump is scattered.

Zoom In Zoom Out Reset image size

If these close encounters, where the minimum distance between the black hole and the clump is ≲4 pc, are the primary source of large fluctuations in the separation of the black holes, then we expect that decreasing ${\epsilon }_{\mathrm{BH}}$ sufficiently would cause these fluctuations to disappear. More precisely, we expect the large fluctuations to disappear once ${F}_{\mathrm{BH}}(\epsilon )/{F}_{\mathrm{cl}}(\epsilon )\gt 1$ for all clumps in the simulation.

In Figures 12–16 we show the evolution of the SMBHs' separation for different values of ${\epsilon }_{\mathrm{BH}}.$ The smallest value that we choose for ${\epsilon }_{\mathrm{BH}}$ is still much greater than the Schwarzschild radius of the black holes in our simulations (${R}_{\mathrm{sch}}\;=\;4.78\times {10}^{-6}\;\mathrm{pc}$). We note that as we change ${\epsilon }_{\mathrm{BH}},$ the orbits of the black holes change. So, in simulations with the same SFR but different ${\epsilon }_{\mathrm{BH}},$ the black holes do not have close encounters with the same gaseous clumps. We found that in three of our simulations (runs C0015_.004, C005_.004, and C05_.004), making ${\epsilon }_{\mathrm{BH}}$ smaller causes the large fluctuations to disappear. As expected, the better resolution of the gravitational force of the black holes in these simulations allows them to disrupt the clumps in all close encounters. In contrast, we found that in the runs C0005_.004 and C015_.004, some large fluctuations in the separation of the pair occurred (see peaks enclosed by red circles in Figures 12 and 15) that are even larger than the fluctuations observed in the simulations with the same SFR, but larger ${\epsilon }_{\mathrm{BH}}$ (runs C0005 and C015). We analyzed each of these large fluctuations to determine why they are present even when the black hole's gravitational force is resolved down to scales of 0.004 pc, where ${F}_{\mathrm{BH}}({\epsilon }_{\mathrm{BH}})$ is larger.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the case of run C0005_.004, we identify two large fluctuations (enclosed by red circles in Figure 12). In Figure 17, we show the evolution of the orbits of the black holes during these two large fluctuations. The upper row of Figure 17 shows the orbits of the black holes during the first of these large fluctuations, and the lower row of Figure 17 shows the second of these large fluctuations. We found that these two large fluctuations occurred because the black holes follow orbits of different radius and eccentricity, but with respect to nearly the same center, which is located close to the densest and most massive gas clump of the simulation (${n}_{\mathrm{cl}}\;=\;7.4\times {10}^{8}\;{\mathrm{cm}}^{-3}$ and ${M}_{\mathrm{cl}}\;=\;8.8\times {10}^{8}\;{M}_{\odot }$). Therefore, the large fluctuations in run C0005_.004 are not the result of close encounters with high-density clumps. Instead they are the result of the movement of the black holes, which, at the moment of these fluctuations, are not bound together.

Zoom In Zoom Out Reset image size

Run C015_.004 also has two large fluctuations in the separation of the black holes (enclosed by red circles in Figure 15). We show the evolution of the orbit of the black holes during these two large fluctuations in Figure 18. From the upper row of Figure 18, we found that the first large fluctuation is the result of scattering of one of the black holes due to gravitational interaction with a high-density gas clump. This gas clump had in turn been ejected from a past scattering with another high-density gas clump. In the lower row of Figure 18, we see that the second large fluctuation is the result of scattering with the same gas clump responsible for the first large fluctuation. Both of these close encounters have a minimum distance that is greater than 10 pc, and therefore they are not artificially produced by low resolution of the gravitational force of the black holes. However, the clump that produces these two large fluctuations in the separation of the black holes is the most dense and massive clump in run C015_.004 ($3.5\times {10}^{8}\;{\mathrm{cm}}^{-3}$ and $4\times {10}^{8}$ ${M}_{\odot }$). So this clump is the most extreme clump in our simulation and hence the two scatterings produced by it are highly unlikely to happen in a real CND. We note that shortly after the last scattering (∼2 Myr), this clump is disrupted by one of the black holes, and after this disruption the orbital decay of the SMBH pair continues relatively smoothly.

Zoom In Zoom Out Reset image size

Intense, fluctuating gravitational torques experienced by SMBH pairs in a clumpy medium have also been observed in numerical simulations by other authors (Escala et al. 2004; Fiacconi et al. 2013; Roskar et al. 2014). To compare our results with the ones obtained by Fiacconi et al. (2013), we estimate the effective density of the black holes and the clumps of their simulations. Using Figure 3 of Fiacconi et al. (2013), we estimate that the density of their clumps ranges between $0.8\times {10}^{6}$ and $13\times {10}^{6}\;{\mathrm{cm}}^{-3}$. This estimate is a lower bound for the density because we assume that the vertical size of their clumps is comparable to the thickness of their disk. If we compare this density range with the effective density of their black holes (which is ${\rho }_{\mathrm{BH}}\;=\;7.$ $9\;{m}_{{\rm{H}}}\;\times {10}^{7}\;{\mathrm{cm}}^{-3}$), we find that ${\rho }_{\mathrm{BH}}$ is typically greater than the density of their clumps and therefore their simulations do not suffer from a badly resolved black holes force. However, from our estimation, we found that the density of these gas clumps is greater or equal than the density of the densest observed molecular clouds. This means that in these simulations the stochastic gravitational torques experienced by the black holes due to the gravitational interaction with the densest gas clumps can be overestimated.

We note that, although in our simulations the maximum density of the gaseous clumps is one or two orders of magnitude greater, the delay on SMBH orbital decay is comparable with that obtained by Fiacconi et al. (2013) for SMBHs that are initially in a circular orbit.

In Fiacconi et al. (2013) they found that the SMBH orbital decay is slower for SMBHs that are initially in a eccentric orbit. This dependence between the orbital decay timescale and the initial eccentricity may be caused by an indirect retarding effect on the SMBH orbital decay that is stronger for SMBHs in an initially eccentric orbit. This indirect effect is caused by the spiral arms produced by the high-density gaseous clumps, which exert gravitational torques on the SMBHs and perturb the hydrodynamical wake of the SMBHs. We expect to find a similar behavior in our simulations; however, as our simulations have star formation included, the clumps' spiral arms may be transformed into stars and their indirect effect on the SMBH orbital decay may be smaller than observed in simulations without star formation. Therefore, it is not clear if in our simulations a SMBH that is initially in an eccentric orbit will have a slower orbital decay than a SMBH that is initially in a circular orbit.