Hot Circumsingle Disks Drive Binary Black Hole Mergers in Active Galactic Nucleus Disks

Active galactic nucleus (AGN) disks have been proposed as one of the most promising locations for producing some of the detected stellar-mass binary black hole (BBH) and neutron star mergers (e.g., McKernan et al. 2012; Bartos et al. 2017; Stone et al. 2017; Leigh et al. 2018; Gröbner et al. 2020; Kimura et al. 2021; Wang et al. 2021; Zhu et al. 2021), especially for the recent detection of the heaviest BBH merger event GW190521 with a possible electromagnetic counterpart (Abbott et al. 2020; Graham et al. 2020). This "AGN channel" could be an intriguing alternative to other traditional formation channels, including isolated binary evolution (e.g., Belczynski et al. 2016; Mandel & de Mink 2016), triple and quadruple systems (e.g., Fernández & Kobayashi 2019; Fragione & Kocsis 2019; Liu & Lai 2021), globular clusters (e.g., Benacquista & Downing 2013; Rodriguez et al. 2016), and nuclear star clusters (e.g., O'Leary et al. 2009; Fragione et al. 2021).

Most studies for binary mergers in AGN disks are based on N-body (e.g., Secunda et al. 2019), Monte Carlo (e.g., Yang et al. 2019; Tagawa et al. 2020; McKernan et al. 2021), semianalytical estimates (Samsing et al. 2020), or hydrodynamical simulations (Baruteau et al. 2011; Kaaz et al. 2021; Li & Lai 2022). Most of these works assume that BBHs contract in AGN disks. However, this is not a certainty, as we have shown in our previous paper (Li et al. 2021, hereafter L21) that prograde BBHs can expand their orbit when the circumsingle disks (CSD) around each BH are adequately resolved. Past studies, such as Baruteau et al. (2011), used excessively large gravitational softenings that effectively removed the important CSD spiral arms. These spirals have been shown to be the key component in driving binary expansion, and so resolving them is crucial (e.g., Tang et al. 2017; Moody et al. 2019; Muñoz et al. 2019; Duffell et al. 2020; Heath & Nixon 2020; Tiede et al. 2020; Dittmann & Ryan 2021; D'Orazio & Duffell 2021; Zrake et al. 2021; L21).

However, L21 assumed a smooth temperature over the whole disk, which may be unrealistic, especially for the regions within the Bondi/Hill radius of each BH. In particular, more realistic treatments of thermodynamics and viscosity are needed to take into account the feedback from the BH. Feedback may be driven by super-Eddington accretion onto stellar-mass BHs. Since luminous AGNs are known to accrete at a few percent of the Eddington rate for supermassive BHs (SMBH) (e.g., Netzer 2015), if an embedded stellar-mass BH accretes at just a fraction of that rate, the resulting accretion rate will be orders of magnitude greater than its Eddington rate.

Numerical MHD simulations for super-Eddington accretion flows show that the midplane disk temperatures around stellar-mass BHs are usually around ∼10⁷ K (e.g., Jiang et al. 2014), which are expected to be much higher than those of the global disk due to higher mass BHs having lower disk temperatures (e.g., Shakura & Sunyaev 1973; Abramowicz et al. 1980).

The heat released by the embedded luminous objects may impact how the binary orbit evolves over time (Szulágyi et al. 2016; Masset 2017; Hankla et al. 2020). Additionally, strong outflows and radiation from accretion disks around BHs could impose significant mechanical and radiative feedback effects on the surrounding environment (Proga & Kallman 2002; Jiang et al. 2014; Yuan et al. 2018). Numerous observational evidence for disk winds has also been accumulated both for cold and hot accretion flows around BHs (e.g., Tombesi et al. 2015; Shi et al. 2021). Those feedback effects could possibly change the thermodynamics around CSDs and thus modify the BBH evolution (Souza Lima et al. 2017; del Valle & Volonteri 2018).

In this work, we extend L21 to study the effect of CSD temperature structure on the BBH orbital evolution in AGN disks using 2D hydrodynamical simulations. To approximate the effect of feedback from each BH, we model the CSD temperature structure with a simple phenomenological model.

This Letter is organized as follows. In Section 2 we introduce the method and models presented in this work. In Section 3 we explore the binary evolution for different CSD temperature profiles. Finally, we summarize the main results of this work and briefly discuss the implications for gravitational wave observations in Section 4.

Most of the model setup and numerical methods are the same as L21. Here we only briefly introduce the key features and defer a more detailed discussion to Appendix A.

We place an equal mass BBH on a prograde, circular orbit around a central SMBH. The orbital radius and orbital frequency of the BBH's center of mass are denoted r₀ and Ω₀, respectively. The BBH is initialized to a circular orbit with semimajor axis a_bin = 0.23R_H, where ${R}_{{\rm{H}}}={\left[{M}_{\mathrm{bin}}/(3{M}_{\mathrm{SMBH}})\right]}^{1/3}{r}_{0}$ is the Hill radius of the BBH, and we set M_bin = 0.002M_SMBH.

The AGN disk is modeled as an isothermal, 2D viscous Shakura & Sunyaev (1973) accretion disk that we initialize to a vertically averaged surface density profile Σ(r) ∝ r^−1/2 and a constant α-viscosity with α = 0.004. We work in a cylindrical, (r, φ), coordinate system centered on the SMBH with radial extent [0.4,2.0] r₀ and uniform resolution equaling ∼110 cells within R_H. The disk feels the gravitational force of the SMBH and BBH, the latter of which is softened with a cubic-spline potential (Kley et al. 2009). Unless otherwise stated, we set the gravitational softening length to _bh = 0.15a_bin. In contrast, the BBH does not feel the gravity of the disk. Our results thus apply to the limit where the disk is low enough in mass that the timescale for the BBH orbit to change is slow compared to the time it takes for the disk to re-equilibrate.

The disk boundary conditions are set to the initial conditions, and we apply wave-killing regions near both boundaries to prevent wave reflections (de Val-Borro et al. 2006). Since we are only interested in the evolution of the BBH separation and not the center of mass motion, the large-scale flow near the boundaries is relatively unimportant. More sophisticated boundary conditions designed to obtain the correct gravitational torque on the center of mass of the BBH (e.g., Dempsey et al. 2020) are thus not required.

To model the effect of the heating from the BBH, we suppose that the gas cooling timescale is fast enough to maintain a temperature profile described by,

$$\begin{eqnarray}&&{ \mathcal R }T({\boldsymbol{r}})={H}^{2}{{\rm{\Omega }}}^{2}\left[1+\sum _{i=1}^{2}f\left(| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{\mathrm{bh},i}| \right)\right],\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal R }$ is the specific gas constant (which we set to unity) and

$$\begin{eqnarray}&&f(\delta r)={T}_{0}{\left(1+\displaystyle \frac{\delta {r}^{2}}{{\epsilon }_{\mathrm{bh}}^{2}}\right)}^{-\tfrac{{\gamma }_{\mathrm{CSD}}}{2}},\end{eqnarray} \tag{ 2 }$$

is a function that specifies the local temperature enhancement around each BH. Far from the BBH, ${ \mathcal R }T(r)\to {H}^{2}{{\rm{\Omega }}}^{2}$, where H = 0.05r₀(r/r₀) is the vertical disk scale height and Ω is the Keplerian frequency around the SMBH. The gas in the vicinity of each BH heats up to an axisymmetric profile following ${ \mathcal R }T(r)\sim (1+{T}_{0}){H}_{0}^{2}{{\rm{\Omega }}}_{0}^{2}/\delta {r}_{i}^{{\gamma }_{\mathrm{CSD}}}$, where δ r_i = ∣ r − r _bh,i ∣. To maintain a finite temperature, we smooth the temperature enhancement to T₀ as the distance to the BH becomes less than the gravitational softening length. We additionally force f(δ r_i ) = 0 when δ r_i ≥ R_H to maintain the background AGN disk temperature on large scales.

In Section 3, we use the hydrodynamical code LA-COMPASS (Li et al. 2005, 2009) to determine how the two-parameter temperature profile of Equation (1) modifies the torque of the disk on the BBH and the corresponding orbital evolution rate ${\dot{a}}_{\mathrm{bin}}$. Previous studies have found that γ_CSD around accreting BHs is usually in the range of ∼0.75−1.0 (Frank et al. 2002; Yuan & Narayan 2014; Jiang et al. 2014). But, an even shallower profile could be attributed to a disk wind (Li et al. 2019; Sun et al. 2019). To cover a wide range of plausible values for γ_CSD, we simulate disks with γ_CSD ranging from 0 to 2.

To keep our simulations relatively simple so as to isolate the effect of the temperature profile, we limit ourselves to studying 2D, locally isothermal, nonaccreting BBHs. Our choice of BBH mass and disk scale height results in a value of R_H = 1.74H. Because of this, 3D effects may be subdominant and our 2D approximation is more realistic than a simulation with R_H < H. Additionally, we neglect accretion onto the BHs to remove an additional parameter from our simulations (the mass removal rate). Finally, our simple temperature model is meant to be a numerically clean approximation to BH feedback. In particular, we are implicitly assuming that any BH feedback is both localized to the BBH Hill sphere and is only in the form of a temperature enhancement around the binary. This simplification is a natural extension of our previous work, and it serves as a first step toward understanding more sophisticated radiation hydrodynamical simulations.

Each simulation is run until ${\dot{a}}_{\mathrm{bin}}$ reaches a quasi-steady state. We discuss the convergence of our simulations in Appendix B (see Figures 4 and 5). Our choice of M_bin, H/r, and α maintains the same gap profile as the more realistic parameters M_bin = 5 × 10⁻⁵ M_SMBH, H/r = 0.01, and α = 0.01 (Lin & Papaloizou 1993; Baruteau et al. 2011; L21).

3.1. Contracting Binaries Have Hot CSDs

Our main conclusion is shown in the first panel of Figure 1 which plots the sign of ${\dot{a}}_{\mathrm{bin}}$ in the parameter space of (T₀, γ_CSD) for all of our simulations. For each γ_CSD, we identify the value of T₀ at which the BBH switches from expansion to contraction. We find that "cold" CSDs result in binary expansion, whereas "hot" CSDs result in binary contraction. The second panel of Figure 1 displays the actual steady-state values of ${\dot{a}}_{\mathrm{bin}}$ for each γ_CSD set as a function of their T₀ values. Noticeably, the magnitude of ${\dot{a}}_{\mathrm{bin}}$ saturates as T₀ increases for each value of γ_CSD. There is a clear boundary separating the hot, contracting binaries from the cold, expanding binaries that lie along the relation ${T}_{0}\approx 2.4\exp ({\gamma }_{\mathrm{CSD}})$ (shown as the gray line in the first panel).

Zoom In Zoom Out Reset image size

In the third panel of Figure 1, we show that this critical T₀ neatly corresponds to a minimum CSD temperature. We plot the radial CSD temperature profiles from Equation (1) for all of our models and color each line blue for expansion and red for contraction. The radial extent of the CSD is ≈0.4a_bin (Artymowicz & Lubow 1994), and we see that at that distance, all of the contracting binaries have a CSD temperature greater than a factor of ∼3 above the background. Thus, when we refer to a "hot" CSD, we are referring to the entire CSD being hotter than the background AGN disk by a factor of ∼3. "Cold" CSDs can still have hot cores, but a steep temperature profile keeps the outer CSD below this threshold and the binary expanding.

Most importantly, we find that even a temperature enhancement as low as T₀ ≳ 5 for a shallow temperature profile (i.e., γ_CSD ≲ 0.75) can alter the BBHs from expansion to contraction. Such a small T₀ could be easily satisfied in a realistic accretion disk environment as discussed in Section 1. This indicates the inspiral of BBHs in AGN disks may be quite common after considering the hot CSD structures. This also suggests the importance of correct thermodynamics of CSDs around the BBH to account for the specific role of AGN disks in BBH merger events.

It should be noted that the softened regions around each BH, which are not modeled accurately, could modify the ${\dot{a}}_{\mathrm{bin}}$ values in Figure 1. We have checked that if we remove the contribution to ${\dot{a}}_{\mathrm{bin}}$ from within 0.5_bh, ⁴ the resulting ${\dot{a}}_{\mathrm{bin}}$ values are approximately the same as those shown in Figure 1. More importantly, the boundary between expansion and contraction remains the same. In fact, the thickness of the gray line in the first panel reflects the change in the boundary when we include or exclude the softening region from the calculation of ${\dot{a}}_{\mathrm{bin}}$.

3.2. Hot CSDs Have Weaker Torques

To understand why a BBH with a hot CSD contracts, we analyze the BH–disk torque ⁵ distribution in this section. We focus on three typical models with γ_CSD = 1 and T₀ values of 0, 10, and 20.

To derive the torque profiles, consider the specific angular momentum of the binary ${{\ell }}_{\mathrm{bin}}={{\boldsymbol{r}}}_{\mathrm{bin}}\times {\dot{{\boldsymbol{r}}}}_{\mathrm{bin}}$, where r _bin = r _bh,1 − r _bh,2 is the BBH separation vector. The binary acceleration is made up of three terms,

$$\begin{eqnarray}&&{\ddot{{\boldsymbol{r}}}}_{\mathrm{bin}}=-\displaystyle \frac{{{GM}}_{\mathrm{bin}}}{| {r}_{b}{| }^{3}}{{\boldsymbol{r}}}_{\mathrm{bin}}+\delta {{\boldsymbol{f}}}_{\mathrm{smbh}}+\delta {{\boldsymbol{f}}}_{\mathrm{disk}}.\end{eqnarray} \tag{ 3 }$$

These are the Keplerian, SMBH, and disk accelerations, respectively. The first two are taken into account in our LA-COMPASS simulations, whereas the disk acceleration, δ f _disk = f _grav,1 − f _grav,2 is measured from the simulation output. A nonzero disk force torques the binary, changing its ℓ_bin at the rate ${\dot{{\ell }}}_{\mathrm{bin}}={{\boldsymbol{r}}}_{\mathrm{bin}}\times \delta {{\boldsymbol{f}}}_{\mathrm{disk}}$. This expression can be rewritten as,

$$\begin{eqnarray}&&{\dot{{\ell }}}_{\mathrm{bin}}=\left(\left({\boldsymbol{r}}-{{\boldsymbol{r}}}_{\mathrm{bh},2}\right)-({\boldsymbol{r}}-{{\boldsymbol{r}}}_{\mathrm{bh},1})\right)\times ({{\boldsymbol{f}}}_{\mathrm{grav},1}-{{\boldsymbol{f}}}_{\mathrm{grav},2}),\end{eqnarray} \tag{ 4 }$$

and because f _grav,i ∝ ( r − r _bh,i ), ${\dot{{\ell }}}_{\mathrm{bin}}$ simplifies to only two terms,

$$\begin{eqnarray}&&{\dot{{\ell }}}_{\mathrm{bin}}=\left({\boldsymbol{r}}-{{\boldsymbol{r}}}_{\mathrm{bh},2}\right)\times {{\boldsymbol{f}}}_{\mathrm{grav},1}+({\boldsymbol{r}}-{{\boldsymbol{r}}}_{\mathrm{bh},1})\times {{\boldsymbol{f}}}_{\mathrm{grav},2}.\end{eqnarray} \tag{ 5 }$$

The first term can be thought of as the torque on BH 1 measured in a coordinate system centered on BH 2, and vice-versa for the second term. In other words, Equation (5) can be written in the more familiar form,

$$\begin{eqnarray}&&{\dot{{\ell }}}_{\mathrm{bin}}=\sum _{j=1}^{2}\int {\rm{\Sigma }}\displaystyle \frac{\partial {{\rm{\Phi }}}_{\mathrm{bh},j}}{\partial {\varphi }_{i}}\delta {r}_{i}d\delta {r}_{i}d{\varphi }_{i}\end{eqnarray} \tag{ 6 }$$

where δ r_i ≡ ∣ r − r _bh,i ∣ and φ_i define a cylindrical coordinate system centered on BH i, and Φ_bh,j is the potential of the companion BH j. As such, we define two torque densities,

$$\begin{eqnarray}&&{{dT}}_{{\rm{g}},1}/{dA}={{\rm{\Sigma }}}_{2}^{{\prime} }\partial {{\rm{\Phi }}}_{\mathrm{bh},1}/\partial {\varphi }_{2}\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}&&{{dT}}_{{\rm{g}},2}/{dA}={{\rm{\Sigma }}}_{1}^{{\prime} }\partial {{\rm{\Phi }}}_{\mathrm{bh},2}/\partial {\varphi }_{1},\end{eqnarray} \tag{ 8 }$$

where ${{\rm{\Sigma }}}_{i}^{{\prime} }$ are the nonaxisymmetric surface densities around each BH and dA is the area element for the appropriate cylindrical coordinate system. The total torque density for the BBH is then ${\dot{{\ell }}}_{\mathrm{bin}}=\int ({{dT}}_{{\rm{g}},1}/{dA}+{{dT}}_{{\rm{g}},2}/{dA}){dA}$.

Figure 2 shows the azimuthally and time-averaged specific torques t_ex,i = δ r_j ∫dT_g,i /dAd φ_j (upper panels), and cumulative torques ${T}_{{\rm{g}},i}={\int }_{0}^{\delta {r}_{j}}{t}_{\mathrm{ex},i}(r){dr}$ (lower panels), acting on BH i centered on BH j. These are time-averaged over 600 BBH orbits once the simulation has reached a steady state. Each column focus on one value of T₀, and for each panel we show both the total (black lines), as well as the individual torques on each BH (green and red lines). The total torque on the binary is the constant value of T_g,1 + T_g,2 at δ r ≫ 3a_bin. This value is positive for T₀ = 0 (indicating expansion), and negative for T₀ = 10 and 20 (indicating contraction). Note that all of the radial profiles are symmetric with respect to each BH. We thus focus only on the torque distribution acting on one of the BHs.

Zoom In Zoom Out Reset image size

We divide the domain into three important regions indicated by the vertical lines. The inner region (IR; δ r ≤ 0.5 a_bin) is the torque on the BH from the companion BH's CSD. This region tends to torque the binary up (t_ex > 0), but is weak overall. Conversely, the outer region (OR; δ r ≥ 1.5 a_bin) tries to torque the binary down (t_ex < 0) and is an important contribution to the overall torque. Lastly, there is the CSD region (0.5 a_bin < δ r < 1.5 a_bin). Gas in this region experiences the strongest torque from the BH, and the total torque from the CSD depends on how asymmetric the torque from the inner CSD region (where t_ex < 0) compares to the torque from the outer CSD region (where t_ex > 0). We see that in most cases, the outer CSD contributes slightly more torque leading to a net positive torque coming from the BH's CSD. The magnitude of the CSD torque, however, becomes weaker overall as the CSD temperature increases. In fact, at T₀ = 10, the positive CSD torque is so weak that the negative OR torque results in binary contracting.

Quantitatively, the torque from the CSD region, normalized to ${{\rm{\Sigma }}}_{0}{r}_{0}^{2}{a}_{\mathrm{bin}}^{2}{{\rm{\Omega }}}_{\mathrm{bin}}^{2}$, goes from 5.0 × 10⁻² when T₀ = 0 to 5.0 × 10⁻³ when T₀ = 10, and even to negative −2.5 × 10⁻³ when T₀ = 20. On the other hand, the torque from the OR region only falls from −4.2 × 10⁻² when T₀ = 0 to −1.8 × 10⁻² when T₀ = 20. As the CSDs get hotter, the magnitude of t_tex decreases overall.

We provide more context for the three regions we have just described in Figure 3 where we show the 2D, time-averaged distributions of surface density and torque density. These profiles are calculated from stacking a large number of snapshots in the rotating frame of the binary (Muñoz et al. 2019; L21), centered on one of the BHs and plotted in Cartesian (x, y) coordinates. From the surface density in the first column and the torque density in the second column we divide the space up into four regions shown by the contours. There are the IR and CSD region defined by the cyan and orange circles, respectively, with radii 0.5a_b centered on each BH. Additionally, there is the new cavity region (denoted CAV) that lies within the red contours. This is meant to approximate the location of the low-density regions between the BHs' CSDs and the outer spiral arms. Finally, everything not enclosed by a contour we place into the OR region. The fourth column shows the cumulative torque profiles for each of these four regions.

Zoom In Zoom Out Reset image size

From the second column, we see that the torque density for each region has a large antisymmetric component about y = 0 and x = x_c , where x_c = 0 for the IR region, x_c = 0.5a_bin for the CAV and OR regions, and x_c = a_bin for the CSD region. The antisymmetric component does not contribute to the total torque (i.e., the green subregions mostly cancel the purple subregions), and so in the third column we show only the symmetric component of the torque for each region. To do this we define the two symmetric torque densities (one about each axis),

$$\begin{eqnarray}\begin{array}{rcl}{\left.\displaystyle \frac{{{dT}}_{g,i}}{{dA}}\right|}^{x} & = & \displaystyle \frac{{{dT}}_{g,i}}{{dA}}(+| x-{x}_{c}| ,y)+\displaystyle \frac{{{dT}}_{g,i}}{{dA}}(-| x-{x}_{c}| ,y),\\ {\left.\displaystyle \frac{{{dT}}_{g,i}}{{dA}}\right|}^{\mathrm{sym}} & = & {\left.\displaystyle \frac{{{dT}}_{g,i}}{{dA}}\right|}^{x}(x,+| y| )+{\left.\displaystyle \frac{{{dT}}_{g,i}}{{dA}}\right|}^{x}(x,-| y| ),\end{array}\end{eqnarray} \tag{ 9 }$$

and then plot the dT_g,i /dA∣^sym for each region in one of the four quadrants in the third column.

Using the symmetric torque and T_g profiles, we come to the following picture of the disk torque on the binary and its dependence on temperature. There is a competition between the positive torques from the IR + CSD regions and the negative torques from the CAV + OR regions. At low temperature enhancements (e.g., T₀ = 0), there is generally more positive torque coming from each BH's CSD that can overcome the CAV and OR torques. When this happens, the relative contribution between the IR and CSD is typically in favor of the CSD. In the T₀ = 0 example, the CSD contributes roughly three times as much positive torque as the IR region.

As T₀ increases, two things happen. First, the spiral arms in the CSDs weaken and their corresponding positive torque decreases. There is roughly a factor of two drop in the IR + CSD torque going from T₀ = 0 to T₀ = 10. Second, the torque from the OR region increases dramatically—almost by a factor of 10 across the same T₀ jump. The torque from the CAV region also increases, but less dramatically (a factor of ∼6 increase). These two effects combined result in a net negative torque on the BBH in the T₀ = 10 case. Going to T₀ = 20 almost completely removes the torque from the BHs' CSDs and only slightly decreases the total torque from the CAV + OR regions.

As a final point, we note that it is not the relative temperatures between the CSDs and circumbinary disks (CBDs) that determine whether the binary expands or contracts. Rather, it is the overall magnitude of the CSD and CBD temperatures that matters. This can be seen by our results at γ_CSD = 0, which raise the CSD and CBD temperatures by the same uniform factor. The fact that these simulations also find that BBHs contract when the CSD is hotter than the background by a factor of ∼3, indicates that the CSD torques weaken enough for the negative CBD torque to set the overall sign of the torque.

We have performed a series of global 2D hydrodynamical simulations of an equal mass, circular BBH embedded in an accretion disk to study the effect of the CSD temperature on the evolution of the binary separation. The CSD temperature profile around each BH is approximated with a phenomenological power-law model designed to mimic the feedback from each BH. We find that when the entire CSD is a factor of ∼3 hotter than the AGN disk, the positive CSD torques are suppressed enough so that the negative outer torques drive the binary to an eventual merger.

A hot CSD can likely arise from strong feedback fueled by super-Eddington accretion onto the BH. This feedback could take the form of jets, winds, and strong radiation (e.g., Jiang et al. 2014; Wang et al. 2021). All of these can modify the flow structures, e.g., whether the CSDs can exist or not, in addition to their effect on the thermodynamics around BBHs. We should point out that we do not consider the accretion onto the BBH in the current work. Although L21 suggests that accretion has a minor effect on the binary orbital evolution for cold CSDs, its effect when coupled with the modified thermodynamics around BBHs still needs to be addressed in the future with detailed simulations.

It is worth noting that almost all simulations of isolated binaries use a temperature profile equal to ${ \mathcal R }T({\boldsymbol{r}})\approx {h}_{0}^{2}{\sum }_{i}{{GM}}_{i}/| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{\mathrm{bh},i}|$, which is similar to Equation (1) with γ_CSD = 1 (e.g., Tang et al. 2017; Moody et al. 2019; Muñoz et al. 2019, 2020; Duffell et al. 2020; Heath & Nixon 2020; Tiede et al. 2020; Dittmann & Ryan 2021; D'Orazio & Duffell 2021; Zrake et al. 2021). The green shaded region in Figure 1 shows the large range of T₀ values used in a sample of those studies. The temperature enhancement in the CSDs can get up to well over a factor of 40 in some cases.

Despite the fact that the CSDs found in isolated binary simulations are "hot," most of these studies still find binary expansion. ⁶ We suspect that this discrepancy lies in the strength of the cavity and outer spiral torques for an isolated binary. In particular, simulations of isolated binaries find much wider and deeper cavities than embedded binaries. This makes the negative cavity and outer region torques on an embedded binary much stronger compared to an isolated binary simply because these regions are "closer" to the binary (i.e., the cavity is shallower and filled in more). This means that in order for an embedded binary to expand, it needs a much stronger CSD torque compared to an isolated binary.

If this picture is correct, the transition point in T₀ for an isolated binary should be at a much larger value than what we find for an embedded binary. This also explains why simulations of isolated binaries in colder disks find contraction (Heath & Nixon 2020; Tiede et al. 2020; Dittmann & Ryan 2022). These studies find that lowering the disk temperature (at fixed viscosity) raises the density near the edge of the cavity, i.e., the peak of the azimuthally averaged surface density profile in the CBD moves closer to the binary (see Figure 10 in Tiede et al. 2020 and Figure 7 in Dittmann & Ryan 2022). Because the cavity wall is closer to the binary and at a higher density, this strengthens the negative torques from this region, making it harder for the binary's CSD torques to push it apart. Indeed, Figure 8 in Tiede et al. (2020) shows that the "cavity" torque becomes more negative as the disk becomes colder.

Our study highlights the important role of thermodynamics around BBHs in their orbital evolution. Because our temperature model is only phenomenological, further work is required to self-consistently determine the CSD temperature profile. In particular, this will require simulations with treatments for radiation and possible feedback from the BH. Calculating the correct feedback mechanism will require simulations near the actual accretion region of the BH which can only be captured by detailed general-relativistic radiative magnetohydrodynamic (GRRMHD) simulations. Such results would then need to be fed into the large-scale simulations presented here—possibly with a subgrid-scale model. The effects of magnetic fields and disk self-gravity may also be important for the structure of the CSDs. Finally, these effects should be included in fully 3D simulations to capture the full BBH–disk interaction problem.

We thank very beneficial comments from Zoltan Haiman and a very constructive report from the referee that helps to clarify some statements of the paper. We gratefully acknowledge the support by LANL/LDRD under project number 20220087DR. This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001. This work is approved for unlimited release under LA-UR-21-31722. Softwares: LA-COMPASS (Li et al. 2005, 2009), Numpy (van der Walt et al. 2011), Scipy (Virtanen et al. 2020), Matplotlib (Hunter 2007).

We briefly describe the equations solved in our hydrodynamical simulations. We numerically solve the continuity and momentum equation for the gaseous disk with an embedded BBH,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}+{\rm{\nabla }}\cdot ({\rm{\Sigma }}{\boldsymbol{v}})=0,\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial ({\rm{\Sigma }}{\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot ({\rm{\Sigma }}{\boldsymbol{v}}{\boldsymbol{v}})=-{\rm{\Sigma }}{\rm{\nabla }}{\rm{\Phi }}-{\rm{\nabla }}P+{\rm{\Sigma }}{{\boldsymbol{f}}}_{\nu },\end{eqnarray} \tag{ A2 }$$

where v is the gas fluid velocity, f _ν is the viscous force from the Shakura & Sunyaev (1973) disk, and P is the gas pressure determined from Equation (1) with an isothermal equation of state. The total gravitational potential of the SMBH-BBH system is

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Phi }} & = & -\displaystyle \frac{{{GM}}_{\mathrm{SMBH}}}{\left|{\boldsymbol{r}}\right|}\\ & & +\displaystyle \sum _{i\,=\,1}^{2}\displaystyle \frac{{M}_{\mathrm{bh},i}}{{M}_{\mathrm{SMBH}}}\left[-\displaystyle \frac{{{GM}}_{\mathrm{SMBH}}}{{D}_{i}}+{{\rm{\Omega }}}_{\mathrm{bh},i}^{2}{{\boldsymbol{r}}}_{\mathrm{bh},i}\cdot {\boldsymbol{r}}\right],\end{array}\end{eqnarray} \tag{ A3 }$$

where Ω_bh,i is the orbital frequency of each BH with respect to the SMBH, and

$$\begin{eqnarray}{D}_{i}^{-1}(\delta {r}_{i})=\left\{\begin{array}{ll}\displaystyle \frac{1}{\delta {r}_{i}}\left[{\left(\displaystyle \frac{\delta {r}_{i}}{{\epsilon }_{\mathrm{bh}}}\right)}^{4}-2{\left(\displaystyle \frac{\delta {r}_{i}}{{\epsilon }_{\mathrm{bh}}}\right)}^{3}+2\displaystyle \frac{\delta {r}_{i}}{{\epsilon }_{\mathrm{bh}}}\right], & \quad \mathrm{if}\ \delta {r}_{i}\leqslant {\epsilon }_{\mathrm{bh}}\\ \displaystyle \frac{1}{\delta {r}_{i}}, & \quad \mathrm{if}\ \delta {r}_{i}\gt {\epsilon }_{\mathrm{bh}}\end{array}\right.\end{eqnarray} \tag{ A4 }$$

determines the softened potential from each BH (Kley et al. 2009). The terms in the bracket of Equation (A3) describe the potential from the BBH in Equation (6).

In this section, we show that our simulations are in steady state by verifying that ${\dot{a}}_{\mathrm{bin}}$ is converged in time, and that the radial mass accretion profile centered on the BBH is spatially constant. The binary semimajor axis evolution rate ${\dot{a}}_{\mathrm{bin}}/{a}_{\mathrm{bin}}$ due to the gravitational force from the gaseous disk is

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{a}}_{\mathrm{bin}}}{{a}_{\mathrm{bin}}}=-\displaystyle \frac{{\dot{\varepsilon }}_{\mathrm{bin}}}{{\varepsilon }_{\mathrm{bin}}},\end{eqnarray} \tag{ B1 }$$

where ε_bin = −GM_bin/(2a_bin) is the specific energy of the binary. We measure the power delivered to the binary by the disk as,

$$\begin{eqnarray}&&{\dot{\varepsilon }}_{\mathrm{bin}}=({\dot{{\boldsymbol{r}}}}_{\mathrm{bh},1}-{\dot{{\boldsymbol{r}}}}_{\mathrm{bh},2})\cdot ({{\boldsymbol{f}}}_{\mathrm{grav},1}-{{\boldsymbol{f}}}_{\mathrm{grav},2}).\end{eqnarray} \tag{ B2 }$$

This is used to compute ${\dot{a}}_{\mathrm{bin}}$ every time-step of the simulation.

The time evolution of ${\dot{a}}_{\mathrm{bin}}/{a}_{\mathrm{bin}}$ for four models is shown in Figure 4. The black line is the model without any enhancement of CSD temperature compared to the global disk (i.e., T₀ = 0), while the blue and red lines correspond to the case of T₀ = 10 and T₀ = 20, respectively. The magenta lines show the model with T₀ = 20 but decreasing the softening scale to 0.075a_bin and increasing the resolution by a factor of two both radially and azimuthally (n_r × n_ϕ = 4096 × 16384). Although the raw data without time-averaging in the upper panel shows strong variabilities, the running window time-averaging plot in the lower panel shows that all ${\dot{a}}_{\mathrm{bin}}/{a}_{\mathrm{bin}}$ values reach a steady state after ∼4000 BBH orbits.

Zoom In Zoom Out Reset image size

As we decrease the softening scale by half, the ${\dot{a}}_{\mathrm{bin}}/{a}_{\mathrm{bin}}$ evolution does not change significantly compared with the two T₀ = 20 models shown in Figure 4, except that the oscillation becomes stronger for the smaller softening case. Our results are thus converged in time, resolution, and softening length.

Figure 5 plots the azimuthally and time-averaged radial profiles of the mass accretion rate in the disk, $\dot{M}(\delta r)\,=2\pi r\left\langle {\rm{\Sigma }}\delta {v}_{r}\right\rangle$, for two γ_CSD = 1 simulations. We compute $\dot{M}$ by stacking the 2D profiles of the radial momentum in the frame centered on the BBH center of mass, Σδv_r , and take the time and azimuthal average. We find that outside of δ r > a_bin, the $\dot{M}$ profiles are both constant and equal to zero. For a nonaccreting BBH, $\dot{M}({R}_{{\rm{H}}})\approx 0$ indicates that the mass in the Hill sphere has reached a steady-state value, and together with Figure 4 verifies that our simulations are converged in time.

Zoom In Zoom Out Reset image size