The Hydrodynamic Evolution of Binary Black Holes Embedded within the Vertically Stratified Disks of Active Galactic Nuclei

Consider a BBH that is embedded within an AGN disk. Assume the center of mass of the binary occupies a circular orbit at distance D from the SMBH, and assume the center of mass corotates with the disk. The disk annuli closer to the SMBH will orbit faster than the BBH, and the annuli farther from the SMBH will orbit slower. In the rest frame of the BBH, this results in a shearing wind that engulfs the binary. Within some sphere of influence, the BBH will feed from the AGN disk, accreting both mass and angular momentum. The first of two relevant length-scales is the Hill radius of the BBH,

$\begin{eqnarray}&&{R}_{{\rm{H}}}=D{\left(\displaystyle \frac{q}{3}\right)}^{1/3},\end{eqnarray} \tag{ 1 }$

where D is the distance from the SMBH to the BBH center of mass, and q = M/M_SMBH is the BBH–SMBH mass ratio. The second length-scale to consider is the Bondi radius,

$\begin{eqnarray}&&{R}_{{\rm{b}}}=\displaystyle \frac{{GM}}{{c}_{{\rm{s}}}^{2}},\end{eqnarray} \tag{ 2 }$

were c_s is the local sound speed at distance D in the disk. At radii larger than R_b, the pressure support of the gas allows it to be unperturbed by the gravity of the BBH. The sphere of influence of the BBH will be limited by the smaller of these two length-scales. If R_H < R_b, the ram pressure of the wind is more important than the gas pressure, and if R_b < R_H, the opposite is true. To parameterize the relative strength of these two regimes, we start by writing down the velocity profile of the wind in the noninertial rest frame of the BBH,

$\begin{eqnarray}&&v(\delta r)={v}_{{\rm{k}}}(D+\delta r)-{v}_{{\rm{k}}}(D),\end{eqnarray} \tag{ 3 }$

where v_k is the Keplerian velocity, and δ r is the relative distance from the BBH. Here, we have essentially assumed that, near (±δ r) the BBH, the velocity profile has a Cartesian geometry. We also neglect the vertical velocity dependence of the flow. It is advantageous to linearize this velocity profile because it will later allow us to model embedded BBHs in a scale-free fashion,

$\begin{eqnarray}\begin{array}{rcl}v(\delta r) & = & \sqrt{\displaystyle \frac{{{GM}}_{\mathrm{SMBH}}}{D+\delta r}}-\sqrt{\displaystyle \frac{{{GM}}_{\mathrm{SMBH}}}{D}}\\ & \approx & \sqrt{\displaystyle \frac{{{GM}}_{\mathrm{SMBH}}}{D}}\left(1-\displaystyle \frac{\delta r}{2D}\right)-\sqrt{\displaystyle \frac{{{GM}}_{\mathrm{SMBH}}}{D}}\\ & \approx & -\displaystyle \frac{1}{2}\delta r{{\rm{\Omega }}}_{{\rm{k}}},\end{array}\end{eqnarray} \tag{ 4 }$

where ${{\rm{\Omega }}}_{{\rm{k}}}=\sqrt{\tfrac{{{GM}}_{\mathrm{SMBH}}}{{D}^{3}}}$ is the Keplerian frequency about the SMBH. The sphere of influence of the BBH is always limited by R_H, and at this radius, the linearized velocity deviates only marginally from the true velocity when q ≪ 1. We define the wind-capture radius, R_w, as the radius at which a wind with the linearized velocity profile v(δ r) is marginally bound to the BBH,

$\begin{eqnarray}&&{R}_{{\rm{w}}}\equiv \displaystyle \frac{2{GM}}{v{({R}_{{\rm{w}}})}^{2}}=2{{Dq}}^{1/3}.\end{eqnarray} \tag{ 5 }$

This quantity differs from the Hill radius only by a constant, i.e., R_w = 24^1/3 R_H. We opt to sometimes use R_w instead of R_H because we are borrowing intuition from Bondi–Hoyle–Lyttleton accretion (e.g., Edgar 2004), where the "accretion radius" describes the region within which a supersonic gas can be gravitationally captured and is defined analagously to our wind-capture radius. Using Equation (4) with δ r = R_w, we define the Mach number of the wind at the wind-capture radius as

$\begin{eqnarray}&&{{ \mathcal M }}_{{\rm{w}}}\equiv \displaystyle \frac{{{\rm{\Omega }}}_{{\rm{k}}}{R}_{{\rm{w}}}}{2{c}_{{\rm{s}}}}.\end{eqnarray} \tag{ 6 }$

Using the relation H/R = c_s/v_k (Pringle 1981), where H/R is the disk aspect ratio, we can rewrite this expression as

$\begin{eqnarray}&&{{ \mathcal M }}_{{\rm{w}}}={q}^{1/3}{(H/R)}^{-1}.\end{eqnarray} \tag{ 7 }$

We have found that our accretion flow is most strongly dictated by ${{ \mathcal M }}_{{\rm{w}}}$. This is illustrated in Figure 1, where we delineate three flow regimes defined by ${{ \mathcal M }}_{{\rm{w}}}$. When ${{ \mathcal M }}_{{\rm{w}}}\ll 1$, then R_b ≪ R_H, and the flow accretes in a quasi-spherical, Bondi-like fashion. When ${{ \mathcal M }}_{{\rm{w}}}\gt 1$, then the wind flowing through the BBH's sphere of influence is supersonic, and has high specific angular momentum. Here, the flow is best characterized as being dominated by the ram pressure of the incident wind. In the intermediate case, such that ${{ \mathcal M }}_{{\rm{w}}}\sim 0.1\mbox{--}1$, the flow still has a large amount of angular momentum, but the higher gas pressure allows the wind to thermalize and produce disk-like structures.

Zoom In Zoom Out Reset image size

It is worth taking a moment to contextualize our wind tunnel framework in the landscape of other accretion flows. Our setup is similar to other local simulations of accretion flows, which have proven effective tools for understanding both single and binary BH accretion in a variety of environments, including the interstellar matter (Antoni et al. 2019; Kaaz et al. 2019, 2021; S. Schrøder et al. 2022, in preparation) and common envelope evolution (MacLeod & Ramirez-Ruiz 2015a, 2015b; MacLeod et al. 2017; Murguia-Berthier et al. 2017; De et al. 2020; Everson et al. 2020). This family of accretion flows essentially considers the permutations of Bondi–Hoyle–Lyttleton accretion (for a review, see Edgar 2004) characterized by the gas properties of the surrounding medium.

The case of the AGN disk is distinguished from these other scenarios by the large angular momentum content available to the embedded binary. In the limiting case of an effectively infinite angular momentum reservoir, one might expect the accretion to resemble the commonly studied thin circumbinary disk ("CBD"; Muñoz et al. 2019, 2020), where BBHs are expected to expand rather than contract due to positive gravitational torques from the surrounding gas. This is similar to the scenario explored by Baruteau et al. (2011) and more recently by Li et al. (2021), who performed global simulations of an AGN disk with an embedded binary system in two dimensions. In Li et al. (2021), they found that prograde embedded binaries also expand, as in the case of CBDs. We emphasize, however, that, for these works to globally resolve the AGN disk and the binary, they are required to use large values of q and model the system in two dimensions. In works such as Li & Lai (2022), who perform two-dimensional shearing box simulations and thus can achieve higher resolutions on the scale of the CBDs and circumsingle disks, they indeed find that the flow behaves differently than in traditional CBDs. In particular, they find that the accretion flow is quasi-steady, exhibiting variability at twice the orbital frequency in the corotating frame. They also find that their binaries contract. However, in the two-dimensional shearing box simulations of Li et al. (2022), they find that the question of binary contraction versus expansion depends on the thermodynamics of the disk, and in the three-dimensional shearing box simulations of Dempsey et al. (2022), they find that that the semimajor axis of the binary is also important.

By using Equation (5), we can write ${R}_{{\rm{w}}}/H\approx 2{q}^{1/3}{(H/R)}^{-1}=2{{ \mathcal M }}_{{\rm{w}}}$. When q is large and H/R is small, R_w ≫ H, which justifies two-dimensional approaches. However, as we will find in the following section, for much of the q − H/R parameter space, H ≳ R_w, and the resulting flow is inherently three-dimensional, supported partially by rotation and partially by pressure gradients. This suggests a dichotomy in the geometry of accretion flows surrounding embedded binaries: in thin accretion disks hosted by lower-mass (≲10⁶–10⁷ M_⊙) SMBHs, the embedded binary may be better represented by the more traditional two-dimensional accretion flows studied by Baruteau et al. (2011) and Li et al. (2021). In either thicker AGN disks (i.e., for high Eddington ratios, or in the outer Toomre-stabilized regions discussed in the next section) or for more massive (≳10⁷ M_⊙) SMBHs, the embedded binaries may be better represented by the vertically stratified accretion flows studied here.

What value of ${{ \mathcal M }}_{{\rm{w}}}$ appropriately represents the flow surrounding an embedded binary? We can answer this by applying our ${{ \mathcal M }}_{{\rm{w}}}$ formalism to physical models of AGN disks. In particular, we consider Shakura–Sunyaev (S-S; Shakura & Sunyaev 1973) and Sirko–Goodman (S-G; Sirko & Goodman 2003) accretion disks. In the classical S-S disk, each disk annulus has Keplerian orbital energy that is gradually dissipated by some unspecified source of viscosity that is parameterized by the quantity α. While the S-S disk forms the basis for most modern disk models, it becomes gravitationally unstable at large distances (typically ≳10²–10³ SMBH gravitational radii). The S-G model addresses these issues by including an unspecified pressure term that forces the disk to be marginally Toomre stable where it would otherwise be unstable. In reality, the disk likely continuously forms stars in this outer region, and the stellar feedback from these populations holds the disk in its marginally stable state.

We can explore how ${{ \mathcal M }}_{{\rm{w}}}$ depends on S-S and S-G AGN disk models by using Equation (7) to relate ${{ \mathcal M }}_{{\rm{w}}}$ to the mass ratio q and disk aspect ratio H/R. The profile of ${{ \mathcal M }}_{{\rm{w}}}$ can depend sensitively on the parameters of the AGN. For simplicity, unless said otherwise, we will assume the host AGN disk has the following parameters for the remainder of this work:

• 1.  
  a 10⁸ M_⊙ SMBH accretes at an Eddington ratio of ${\lambda }_{\mathrm{Edd}}^{(\mathrm{SMBH})}=0.5$, with a radiative efficiency of η = 0.1,
• 2.  
  viscosity is proportional to the total (gas + radiation) pressure, and α = 0.01,
• 3.  
  the embedded, equal-mass binary has a total mass M = 60 M_⊙.

In the case of the S-G model, this is the same disk model used to study the orbital migration of embedded BHs in Secunda et al. (2019), who found that their migrating BHs became trapped at roughly ≈10⁻³–10⁻² pc.

In the left panel of Figure 2, we plot ${{ \mathcal M }}_{{\rm{w}}}$ as a function of distance D in the AGN disk for both profiles. For the S-S profile, we show both a disk with H/R set by gas pressure (solid line) and by radiation pressure (dashed line; Frank et al. 2002). We also shade two regions:

• 1.  
  The lower orange region we designate as Bondi-like (${{ \mathcal M }}_{{\rm{w}}}\lesssim 0.1$). Within this region, $\dot{M}$ (Equation (12)) deviates from the Bondi accretion rate (Bondi 1952) by only ≲1%. If this is the case, the velocity shear in the wind is mostly negligible, and the flow geometry is quasi-spherical.
• 2.  
  The upper blue region is where our binary opens a gap within the AGN disk. This is defined by the gap-opening criterion provided in Equation (15) of Crida et al. (2006), and while it is insensitive to the specific disk profile it does depend on α.

We notice a few features immediately: First, it is unusual for ${{ \mathcal M }}_{{\rm{w}}}$ to exceed unity. An exception is the gas pressure dominated S-S curve, but only at small radii where the disk is in fact radiation pressure dominated at this value of ${\lambda }_{\mathrm{Edd}}^{(\mathrm{SMBH})}$. The radiation pressure dominated S-S curve has large ${{ \mathcal M }}_{{\rm{w}}}$ at large radii, but here the flow is gas pressure dominated. For the S-G profile, which is more credible in the outer regions where the S-S disk is Toomre unstable, the flow is actually Bondi-like at regions beyond ≈1 pc. We note the presence of the inflection point at ≈0.05 pc in the S-G profile, which marks the transition between the inner disk and the Toomre-stabilized outer disk (see for reference Figure 2 of Sirko & Goodman 2003, which depicts the same profiles used here).

Zoom In Zoom Out Reset image size

In the right panel of Figure 2, we plot contours of ${{ \mathcal M }}_{{\rm{w}}}$ in the q–H/R plane (Equation (7)), and shade the same Bondi-like (orange) and gap-opening (blue) regions as in the left panel. The vertical axis in this panel is alternatively labeled in terms of M_SMBH, given M = 60 M_⊙. We have labeled ${{ \mathcal M }}_{{\rm{w}}}=0.5,1$ and 2 isocontours, which correspond to the hydrodynamic simulations that we present in Section 3. The gap-opening criterion is satisfied only for particularly thin accretion disks with particularly high mass ratios. This point deserves emphasis, as there are conflicting claims in the literature regarding gap formation in AGN disks. In Baruteau et al. (2011), the gap formation happens primarily because of their choice of α = 8 × 10⁻⁴ (α ≈ 0.01–0.3 is more realistic for AGN disks, e.g., King et al. 2007), and the gap opening is clearly sensitive to this parameter as can be seen in Figure 2. In Bartos et al. 2017, they also invoke gap opening, but use a disk model that is thinner than those from both S-G and S-S models, causing ${{ \mathcal M }}_{{\rm{w}}}$ to be larger. In Stone et al. (2017), they also find that gaps are unlikely to open anywhere in the disk, and they use the AGN disk model from Thompson et al. (2005). This model is comparable to the S-G model because they also construct profiles that are marginally Toomre stable in the outer regions. Still, how significantly gap-opening affects the orbital evolution of embedded binaries remains unclear. For instance, in Li et al. (2021), they find that the gap depth weakly affects the evolution of the semimajor axis, particularly at larger H/R. This is because the binary evolution is mainly dictated by gas within the Hill sphere(s), which can be readily replenished even when a gap is opened.

To test our ${{ \mathcal M }}_{{\rm{w}}}$ formalism, we performed three-dimensional hydrodynamic simulations of accretion onto embedded BBHs using version 4.5 of the grid-based, adaptive mesh refinement hydrodynamics code FLASH (Fryxell et al. 2000). We use a local, wind tunnel computational domain that is continuously replenished with a shearing gaseous wind (see the arrows in the inset panels of Figure 1). Physically, the velocity profile of the wind is an approximation derived from the Keplerian velocity about the SMBH in the noninertial rest frame of the BBH, which we acquire using the linearized velocity profile defined in Equation (4). For simplicity, we neglect gradients in the gas density and sound speed of the wind, setting them to ambient values ρ_∞ and c_∞ respectively. We evolve our simulations in a scale-free fashion, where c_∞ = R_b = 1, and ρ_∞ = 10. We relate the linearized velocity profile in Equation (4) to our scale-free setup by rewriting it as

$\begin{eqnarray}&&v(y)=y\displaystyle \frac{{c}_{\infty }{{ \mathcal M }}_{{\rm{w}}}^{3}}{2{R}_{{\rm{b}}}}\end{eqnarray} \tag{ 8 }$

where we have substituted δ r for y.

In this fiducial work, we solve the equations of inviscid hydrodynamics rather than the shearing box equations, and leave the latter for a later work. This means that we neglect centrifugal and coriolis accelerations from the SMBH potential. Our initial conditions in the absence of the BBH are setup to be in a steady state without these forces. When constructing the domain, we also neglect the vertical gradients in the wind, acknowledging that when ${{ \mathcal M }}_{{\rm{w}}}\gt 1$ the AGN disk scale height H becomes comparable to the BBH sphere of influence. We also neglect the breaking of the symmetric velocity profile by pressure gradients in the AGN disk. These pressure gradients however modify the velocity profile on the order of (H/R)² (Pringle 1981), which is negligible for much of the parameter space (e.g., Figure 2).

These neglected effects are negligible when a ≪ R_H (where a is the semimajor axis of the binary), but would affect our results when a ≳ R_H. In a real system, binaries with a > R_H are unbound. Some of our simulation parameter choices correspond to this regime, but remain valuable for establishing hydrodynamic scaling relations within our wind tunnel setup. We also neglect the self-gravity of the gas; however, even if self-gravity is important for the equilibrium state of the Toomre-unstable regions of an AGN disk, within our domain the gravity of the binary dominates, and the self-gravity of the gas is dynamically unimportant. The width of our computational domain is 32 R_b in each Cartesian direction, with the wind oriented in the $\pm \hat{x}$ direction and generated in every time step at the x boundaries of the domain. For the y- and z-directions, we use outflow boundary conditions. Our simulation domain size is much larger than the BBH sphere of influence, and consequentially the overall gas mass in the simulation domain remains roughly constant.

It is well known that hierarchical triple systems are sometimes unstable (Eggleton & Kiseleva 1995; Mardling & Aarseth 2001; He & Petrovich 2018). The results of Eggleton & Kiseleva (1995) suggest that, for circular orbits, the stability of the inner binary is dictated by the semimajor axes of the inner (a) and outer (D) orbits and the mass ratio (q) of the inner binary to the outer massive object (the SMBH, in our case). So, while the hydrodynamics of our simulations are self-similar, they are not necessarily dynamically stable for all parameter choices. We acknowledge this as a limitation of our local simulations, but emphasize that the purpose of this work is to isolate the hydrodynamics of accreting, embedded BBHs. This allows us to develop scaling relations that can then be applied to the dynamically stable parameter space.

We evolve the BBHs using the active sink particles module that was developed for FLASH by Federrath et al. (2010), and use the same implementation of it as Antoni et al. (2019). While we refer the reader to these works for complete descriptions of our sink prescription, the most salient details are as follows. The binary is initialized with its center of mass at x = y = z = 0 and occupies a circular orbit with semimajor axis a in the the x–y plane. The +z-axis is parallel to the angular momentum vectors of the binary orbit and the AGN disk. We model the BBHs using Newtonian gravitational potentials with an absorbing boundary condition centered at the location of each binary companion. These absorbing boundaries have sink radii R_s = 0.0125 R_b, and any gas that crosses this boundary is removed from the domain, with the mass and angular momentum of the removed gas being recorded and added to the corresponding binary member. We also check the convergence of our mass accretion and inspiral rates with the sink radius in Appendix A. While the rate of inspiral is converged, the mass accretion rate decreases with decreasing sink radius. The nonconvergence of the mass accretion rate with the sink radius is a known issue with the sink prescription in Bondi-like accretion flows (Xu & Stone 2019; De et al. 2020), and so the accretion rates derived from these simulations should be taken as upper limits. We adaptively refined around the region closest to the binary members, with smallest cell size $\delta {r}_{\min }\approx 0.0039\,{R}_{{\rm{b}}}$. This means that there are ${R}_{{\rm{s}}}/\delta {r}_{\min }=3.2$ cells across a sink radius, a value that was tested for convergence in Antoni et al. (2019).

In total, we present six simulations, varying ${{ \mathcal M }}_{{\rm{w}}}$ and a/R_b, where a is the BBH semimajor axis. We use values ${{ \mathcal M }}_{{\rm{w}}}=0.5,1,$ and 2, and a/R_b = 0.1 and 1. We use a gamma-law equation of state, and assuming that the gas in the AGN disk can cool efficiently, take γ = 1.1. ⁶ In each simulation, the BBH occupies a circular, equatorial orbit that is prograde with the angular momentum of the wind. Each simulation was run until t = 50 R_b/c_∞, well into the steady state as shown in Appendix A. In the following two subsections, we present the results of these simulations.

In Figure 3, we show the large-scale flow morphology of each simulation. In the top 2 × 3 panel subfigure, we depict 12 × 12 R_b equatorial slices of the gas density. When ${{ \mathcal M }}_{{\rm{w}}}=0.5$, the flow is structured and laminar, with an extended density profile that is wound in a direction prograde with the binary angular momentum. When ${{ \mathcal M }}_{{\rm{w}}}=2$, the flow becomes unwound, with tails that follow the direction of the wind. When ${{ \mathcal M }}_{{\rm{w}}}=1$, the flow shares characteristics of both the ${{ \mathcal M }}_{{\rm{w}}}=0.5$ and 2 simulations. It is worth taking a moment to compare these structures to those found in the global simulations of Li et al. (2021), Baruteau et al. (2011). In these works, the disk is cool and two-dimensional, and spirals around the mini-disks of the binaries are formed and become unwound as they transition into the global AGN disk. While our simulations share some of these features, particularly in the bottom left and middle panels of Figure 3, we generally find much less pronounced spiral patterns. We hypothesize that this is due to the following distinguishing features in our accretion flow. First, the disks in our simulations are hot and partially pressure-supported, making the angular velocity profile sub-Keplerian. This, in turn, will change the radial spacing of the Lindblad resonances, which is responsible for the formation of spiral density waves. Second, as seen most prominently at ${{ \mathcal M }}_{{\rm{w}}}\geqslant 1$, our flow has a clumpy, turbulent density distribution, which will damp any spiral density waves that are formed. Disk-like structures are evident in our two-dimensional slices, and can be quantified with the circularization radius of the wind,

$\begin{eqnarray}&&{R}_{\mathrm{circ}}=\displaystyle \frac{{GM}}{{v}^{2}(y={R}_{\mathrm{circ}})}\to {R}_{\mathrm{circ}}={{ \mathcal M }}_{{\rm{w}}}^{-2}{4}^{1/3}{R}_{{\rm{b}}}.\end{eqnarray} \tag{ 9 }$

For ${{ \mathcal M }}_{{\rm{w}}}=0.5$, 1, and 2, the corresponding circularization radii are R_circ ≈ 6.4, 1.6, and 0.4 R_b. These values roughly correspond to the extent of the density structures in each simulation, particularly when a/R_b = 0.1.

Zoom In Zoom Out Reset image size

The main difference between the a/R_b = 0.1 and 1 simulations is that, in the former, the accretion flow forms a high-density envelope shared by both binary companions, while in the latter each companion has its own envelope. This makes sense, because, when a/R_b = 0.1, the semimajor axis of the binary is much smaller than the sphere of influence, while, when a/R_b = 1, they are comparable.

In the bottom 1 × 3 panel subfigure of Figure 3, we show the three-dimensional gas density contours from each a/R_b = 1 simulation. As also seen in the two-dimensional subfigure, when a/R_b = 1, each binary companion has its own high-density envelope, which overlap more when ${{ \mathcal M }}_{{\rm{w}}}$ is smaller. As ${{ \mathcal M }}_{{\rm{w}}}$ increases, the peak density remains similar, but each density envelope becomes truncated. While in the two-dimensional subfigure, the flow appears disk-like; we see from the three-dimensional panel that the density envelopes are as extended in the vertical directions as they are in the equatorial directions. This suggests that the flow structures surrounding embedded BBHs are partially Bondi-like, and partially disk-like.

The orbits of our simulated BBHs evolve due to a combination of gas drag and the accretion of mass and angular momentum. To compare to our numerical results, we will begin by writing down the analytic binary inspiral rate, $\dot{a}/a$. The full expression for $\dot{a}/a$ is

$\begin{eqnarray}&&\displaystyle \frac{\dot{a}}{a}=2{\left(\displaystyle \frac{\dot{L}}{L}\right)}_{\mathrm{accretion}}-\displaystyle \frac{3\dot{M}}{M}+{\left(\displaystyle \frac{\dot{a}}{a}\right)}_{\mathrm{drag}},\end{eqnarray} \tag{ 10 }$

which assumes that the specific angular momentum of the binary is Keplerian $\left({\ell }=\sqrt{{GMa}}\right)$. We have split the angular momentum evolution into two terms for accretion and gas drag, i.e., $2\dot{L}/L=2{\left(\dot{L}/L\right)}_{\mathrm{accretion}}+2{\left(\dot{L}/L\right)}_{\mathrm{drag}}$, where we have written $2{\left(\dot{L}/L\right)}_{\mathrm{drag}}$ as ${\left(\tfrac{\dot{a}}{a}\right)}_{\mathrm{drag}}$ to be consistent with our analytic estimate of drag, which we will introduce shortly. The contribution from gas drag is found to always be negative in our work. This is in contrast to hydrodynamic simulations of thin, gap-forming circumbinary disks where the drag term is often found to be positive, e.g., Muñoz et al. (2019). We are mainly interested in the latter two terms for mass accretion rate and drag. Note, however, we do include a time series for Equation (10) in Figure 8. To determine $\dot{M}$ analytically, we integrate over streamlines entering the gravitational sphere of influence of the BBH. We start by defining the length-scale,

$\begin{eqnarray}&&{R}_{\mathrm{bw}}=\displaystyle \frac{{R}_{{\rm{b}}}}{1+\tfrac{1}{2}{{ \mathcal M }}_{{\rm{w}}}^{2}},\end{eqnarray} \tag{ 11 }$

which smoothly transitions between R_b when ${{ \mathcal M }}_{{\rm{w}}}\ll 1$, and R_w when ${{ \mathcal M }}_{{\rm{w}}}\gg 1$. This allows us to continuously characterize the sphere of influence of the BBH across flow regimes analytically. The integral for the mass accretion rate is

$\begin{eqnarray}\begin{array}{rcl}\dot{M} & = & f\times {\rho }_{\infty }{\displaystyle \int }_{r\lt {R}_{\mathrm{bw}}}v(y){dA}\\ & = & f\times 4\pi {R}_{\mathrm{bw}}^{2}{\rho }_{\infty }{c}_{\infty }\left(\displaystyle \frac{{{ \mathcal M }}_{{\rm{w}}}^{3}}{8+4{{ \mathcal M }}_{{\rm{w}}}^{2}}+1\right),\end{array}\end{eqnarray} \tag{ 12 }$

where we have included an ad hoc numerical prefactor f. We do this because these calculations assume that all streamlines entering the BBH sphere of influence are accreted; though in reality some fraction of this material is advected away due to gas and ram pressure, such that f < 1. To estimate ${\left(\tfrac{\dot{a}}{a}\right)}_{\mathrm{drag}}$, we follow the approach of Antoni et al. (2019), who determined the inspiral rate for a BBH embedded in a Bondi–Hoyle wind tunnel (e.g., Edgar 2004; Blondin & Raymer 2012). In Appendix B, we modify their derivation for our wind tunnel setup and derive the expression

$\begin{eqnarray}&&{\left(\displaystyle \frac{\dot{a}}{a}\right)}_{\mathrm{drag}}=s\times \displaystyle \frac{8\pi {Ga}{\rho }_{\infty }}{{({GM}/4a+{c}_{\infty }^{2}+{{ \mathcal M }}_{{\rm{w}}}^{2}{c}_{\infty }^{2})}^{1/2}}{\left(1+\displaystyle \frac{2{R}_{\mathrm{bw}}}{a}\right)}^{p}.\end{eqnarray} \tag{ 13 }$

Here, s and p are numerically determined parameters that we will fit for. This calculation assumes that the orbital energy of the binary is dissipated by a gravitational wake induced by each component of the binary as they orbit through a dense gaseous envelope. These wake-induced gravitational torques are different than those that dictate the orbital angular momentum evolution in thin circumbinary disks, which is mainly driven by Lindblad torques from the inner and outer disks. Our main argument for neglecting the study of Lindblad torques is that our accretion flows are pressure dominated, rather than rotation dominated, and we do not have any circumbinary disks. It is possible that there is some transition that occurs at higher ${{ \mathcal M }}_{{\rm{w}}}$, where the angular momentum evolution of the binary transitions to a structure that better resembles the usual circumbinary disk.

In Figure 4, we assess how well our analytic estimates compare to our simulated values for the accretion and drag rates as a function of ${{ \mathcal M }}_{{\rm{w}}}$. Our simulated values of mass accretion rate are calculated by taking the time derivative of the mass of the binary, which is updated each time step during our simulation. The drag rate is determined by recording the net drag force that the gas exerts on the binary at each time step. We express our rates in units of ${t}_{\mathrm{Bondi}}^{-1}={c}_{\infty }^{3}/4\pi {G}^{2}M{\rho }_{\infty }$, which is the Bondi accretion rate per unit mass (Bondi 1952). We have numerically determined that f = 0.58 (Equation (12)), s = 0.35, and p = −0.86 (Equation (13)). Using these adjustments, Figure 4 shows that our analytic estimates for our mass accretion and drag rates reasonably reproduce our simulated results. In general, the binaries accrete and inspiral on similar timescales, with inspiral occurring faster (slower) at larger (smaller) separations. This is in contrast to what is seen in Bondi–Hoyle accretion onto binaries, e.g., Figure 11 of Antoni et al. (2019), where inspiral happens about 3–4 times faster than accretion. We note that, when ${{ \mathcal M }}_{{\rm{w}}}\gg 1$, our local wind tunnel approximation breaks down because the BBH can begin influencing the AGN disk, so care should be taken in extrapolating to larger values of ${{ \mathcal M }}_{{\rm{w}}}$. In those cases, the flow geometry is likely better represented by the results of Baruteau et al. (2011). In contrast, extrapolating our mass accretion prescription to ${{ \mathcal M }}_{{\rm{w}}}=0$ provides no issue, as the flow approaches a Bondi solution in that limit.

Zoom In Zoom Out Reset image size

To determine the evolution of embedded BBHs, it is important to understand how they are first captured. The BBH is either born in the outer regions of the disk or, having previously existed in the NC when the disk was formed, it gradually aligns with the disk. The initial position of the embedded BBH is determined by how it became embedded, which strongly impacts the binary's subsequent evolution.

Roughly ≈80% of stars in the NC are expected to have formed in situ (Antonini et al. 2015). If the BBH is formed directly in the disk (as studied in Stone et al. 2017), it by definition will start in the Toomre-stabilized region as this is where star formation occurs. In Figure 5 of Thompson et al. (2005), they generally find that the star formation rate in nuclear starburst disks increases as a function of radius. This would suggest that the vast majority of embedded BBHs formed in situ will begin their journeys at large (≳1–10 pc) radii in the AGN disk. The possibility that they merge in the inner regions of the AGN disk depends on the relative efficiency of inspiral versus migration.

The dynamical capture of BBHs by AGN disks depends on the distribution of stellar-mass BHs (sBHs) in the NC. This distribution is often chosen by invoking the inferred distribution of the BH cusp surrounding Sag A*, which roughly scales as ∝D^−2.5 (Alexander & Hopman 2009; Bartko et al. 2009) and likely has ∼1000 sBHs in the inner 0.1 pc (Antonini 2014; Hailey et al. 2018). However, the formation of this cusp requires that the stellar population has relaxed around the central SMBH (Bahcall & Wolf 1976). The relaxation timescale for a 10⁶ M_⊙ SMBH can range from 0.1–10 Gyr (O'Leary et al. 2009), and can exceed a Hubble time for more massive SMBHs. So, it may be that heavier SMBHs harbor no BH cusp during the AGN phase, which makes the viability of the dynamical capture channel less certain. For this reason, we favor the in situ formation of embedded BBHs, which should begin their journey in the AGN disk at larger (≳1 pc) distances.

Once BBHs become embedded in the AGN disk, they will simultaneously inspiral, migrate, and grow in mass. In this section, we model these processes by evolving a system of coupled ordinary differential equations for $\dot{D}$ (migration), $\dot{M}$ (accretion), and $\dot{a}$ (inspiral). To guide our intuition, we begin by estimating the timescales associated with each of these processes. We use an S-G AGN disk model with M_SMBH = 10⁸, α = 0.01, ${\lambda }_{\mathrm{Edd}}^{(\mathrm{SMBH})}=0.5$and assume that the viscosity is proportional to the total (gas + radiation) pressure. At a characteristic distance of 1 pc for a 60 M_⊙ binary, the migration timescale is

$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{migr}} & \approx & 61\,\mathrm{Myr}{\left(\displaystyle \frac{0.06}{H/R}\right)}^{-2}{\left(\displaystyle \frac{60\,{M}_{\odot }}{{M}_{\mathrm{BBH}}}\right)}^{-1}{\left(\displaystyle \frac{{10}^{8}\,{M}_{\odot }}{{M}_{\mathrm{SMBH}}}\right)}^{3/2}\\ & & \times {\left(\displaystyle \frac{0.5\,\mathrm{pc}}{D}\right)}^{-1/2}{\left(\displaystyle \frac{{\rm{\Sigma }}}{1540\,{\rm{g}}\,{\mathrm{cm}}^{-2}}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 14 }$

which we acquire from Equation (10) of Paardekooper (2014). This characteristic timescale is the comparable lifetime of the host disk (≈50–150 Myr), suggesting that embedded binaries may migrate significantly (e,g, Secunda et al. 2019, 2020). This equation is for the static torque usually associated with Type I migration, which is valid under our assumption that the embedded BBH only negligibly affects the evolution of the host AGN disk. This expression also neglects self-gravity of the host disk and radiative torques, which may enhance the migration rate (Kley & Nelson 2012).

At large semimajor axes, before gravitational-wave emission dominates, the binaries inspiral due to the combined influence of mass accretion and drag. At distances of ≳0.5 pc, the accretion flow is Bondi-like (e.g., Figure 2), for which the mass doubling time ($\equiv M/{\dot{M}}_{\mathrm{Bondi}}$) is characteristically

$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{Bondi}} & \approx & 148\,\mathrm{Kyr}{\left(\displaystyle \frac{M}{60\,{M}_{\odot }}\right)}^{-1}\left(\displaystyle \frac{\rho }{8\times {10}^{-15}\,{\rm{g}}\,{\mathrm{cm}}^{-1}}\right)\\ & & \times {\left(\displaystyle \frac{{c}_{s}}{54\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)}^{-3}.\end{array}\end{eqnarray} \tag{ 15 }$

This is extremely rapid, which at first glance suggests the embedded binaries should easily balloon in mass. The associated drag timescale (≡the inverse of Equation (13), where for simplicity we take ${{ \mathcal M }}_{{\rm{w}}}=0$, and a = 100 au) is ≈18 Kyr, increasing to ≈72 Kyr at a = 10 au. Together, these three timescales form a hierarchy: the binary inspirals first (mainly due to drag, but also due to mass accretion), then increases in mass, then migrates.

However, the mass accretion rate associated with Equation (15) is super-Eddington throughout much of the disk. This can be seen in Figure 5, where we plot the Eddington ratio of the embedded binary, ${\lambda }_{\mathrm{Edd}}^{(\mathrm{BBH})}$, as a function of distance in the host AGN disk. This figure shows that the embedded binary can be supplied gas at rates of up to ∼10⁴–10⁵ the Eddington ratio. These profiles peak at radii of ∼0.05–0.07 pc, where there is a transition to the star-forming outer regions of the host AGN disk. In practice, we limit the mass accretion rate to its Eddington-limited value,

$\begin{eqnarray}&&{\dot{M}}_{\mathrm{Edd}}=\displaystyle \frac{{L}_{\mathrm{Edd}}}{\eta {c}^{2}}.\end{eqnarray} \tag{ 16 }$

Simulations of super-Eddington accretion have shown that the mass accretion rates can significantly exceed the Eddington rate (McKinney et al. 2014, 2015; Dai et al. 2018), which is essentially set by how small the radiative efficiency (η) is. Thin disks have typical values of η = 0.1, which lowers to roughly η = 0.01 near the Eddington limit, and can be even smaller at significantly super-Eddington rates. For simplicity, we adopt a constant value of η = 0.01. The associated Eddington accretion timescale ($\equiv M/{\dot{M}}_{\mathrm{Edd}}$) is independent of mass,

$\begin{eqnarray}&&{t}_{\mathrm{Edd}}\approx 4.5\,\mathrm{Myr}\left(\displaystyle \frac{\eta }{0.01}\right).\end{eqnarray} \tag{ 17 }$

Zoom In Zoom Out Reset image size

Even with Eddington-limiting our mass accretion rate, it still occurs on a shorter timescale than that in migration. We must make a few modifications for the inspiral rate in our binary evolution prescriptions, as well. First, we include a term for the gravitational-wave inspiral rate of an equal-mass binary (Peters 1964). Second, we neglect the ${(\dot{L}/L)}_{\mathrm{accretion}}$ term in Equation (10), because in real systems this term can only contribute marginally (${ \mathcal O }(\dot{M}{l}_{\mathrm{isco}})$, where l_isco is the specific angular momentum at the innermost stable circular orbit) to the binary's evolution. We also make the assumption that whenever the mass accretion rate is Eddington limited, so is the drag rate, because an Eddington-limited disk will be depleted of the gas that drives the binary to inspiral.

We evolve our binaries in Figure 6, where we show how the binary separation, mass, and AGN disk position change in time. We initialize our evolved binaries with masses M = m₁ + m₂ = 20 (dashed lines) and 60 M_⊙ (solid lines), and initial distances D = 0.5, 1, and 2 pc. We initialize our binary semimajor axes such that they are marginally stable to the hierarchical triple instability outlined in Eggleton & Kiseleva (1995). We use ${\lambda }_{\mathrm{Edd}}^{(\mathrm{SMBH})}=0.5$ on the left subfigure, and ${\lambda }_{\mathrm{Edd}}^{(\mathrm{SMBH})}=0.05$ on the right. Each integration is terminated when the semimajor axis reaches a = 0.1 au, after which the inspiral becomes rapid due to gravitational-wave decay. In all cases, the full evolution occurs on timescales of ≲1–30 Myr. The relative hierarchy of timescales follows our expectations from our analytic predictions; the binary inspirals the quickest, and is driven to inspiral by a combination of drag, gravitational-wave radiation, and mass accretion. Second, the masses begin to increase and, in some cases, double. For each profile, the binaries migrate only marginally before they merge. To highlight uncertainties in our binary evolution prescriptions, we reran each profile with the prescribed mass accretion rate reduced or increased by a factor of 10, the results of which are shown by the error bars in our middle panel.

Zoom In Zoom Out Reset image size

In general, the evolution of the binary happens well within the lifetime of the AGN disk (≈50–150 Myr, labeled τ_AGN in Figure 6). However, a potentially important caveat in this evolution is that we keep the AGN disk static. In reality, the AGN disk can "flicker" on short timescales (≈0.1–10 Myr, labeled δ τ_AGN in Figure 6; e.g., Schawinski et al. 2015). A similar phenomena is seen in the simulations of Angles-Alcazar et al. (2020), who found that the gas supplying luminous quasars can restructure and change the orientation on timescales of ≈0.1–1 Myr. If the binary becomes embedded, but then within δ τ_AGN, the disk restructures and reorients itself; the BBH may once again require its inclination to be damped with respect to the restructured disk for a gas-assisted inspiral to proceed.

The results of Figure 6 have profound implications for the gravitational-wave signals observed by LIGO and Virgo. If a BBH becomes embedded within an AGN disk, it can be driven to merger on timescales ≲1–30 Myr. While the rates at which BBHs form within AGN disks are unclear, our calculations suggest that a gas-assisted inspiral is a viable formation channel for BBH merger events. Furthermore, the BHs in these binaries accrete at significant rates, suggesting that a pair of 30 + 30 M_⊙ BHs can potentially evolve into the pair-instability supernova mass gap. Such a system was recently observed by LIGO/Virgo (Abbott et al. 2020, 2020), with a reported possible association with an EM transient event in an AGN (Graham et al. 2020). If further analysis confirms the association or if additional associations are identified with future LIGO events, such signals would amount to a confirmation of the merger channel studied here.

A caveat is that the accreting BHs in our paradigm should significantly spin up and become aligned, leading to a large effective spin parameter, χ_eff ≳ 0.9. In contrast, the GW190521 gravitational-wave signal was associated with a much smaller value of χ_eff = 0.08. This suggests that the BHs either had low spins individually or had large spins that were misaligned with the orbit. If GW190521 did come from the AGN disk formation channel, then there must be some difference in the evolutionary process that we characterized in Section 4.2. A few plausible scenarios are as follows:

• 1.  
  While we focus on embedded BBHs, individual BHs will migrate and accrete mass from the disk in an analogous way. The inner parsec surrounding the host SMBH may be replete with these single BHs, each of which can potentially grow to sufficiently large masses that they enter the PISN mass gap. These BHs could then pair and merge dynamically later in their life. In this scenario, the BH spins would be randomly oriented, potentially leading to a low χ_eff, regardless of the individual BH spins.
• 2.  
  The scenario presented here involves no dynamical interactions while the BBH evolves within the disk. We can calculate the interaction time between a BBH and a tertiary star with a randomly oriented orbit:

  $\begin{eqnarray*}\begin{array}{rcl}t & = & \displaystyle \frac{1}{n\sigma V}\\ & \simeq & 2{\left(\displaystyle \frac{n}{{10}^{3}\ {\mathrm{pc}}^{-3}}\right)}^{-1}{\left(\displaystyle \frac{a}{100\ \mathrm{au}}\right)}^{-2}{\left(\displaystyle \frac{{v}_{{\rm{k}}}}{650\ \mathrm{km}\ {{\rm{s}}}^{-1}}\right)}^{-1}\mathrm{Myr},\end{array}\end{eqnarray*}$

  where we have adopted a cross section, σ, derived from the orbital separation, and an interaction velocity V corresponding to the Keplerian orbit of a BBH at a distance of ∼1 pc from a 10⁸ M_⊙ SMBH. For a stellar density of 10³ pc⁻³ and an orbital separation of 10² au, we find an interaction time of ≈2 Myr, shorter than a merger time of ≈10 Myr. These strong dynamical encounters can fundamentally alter the BBHs orbit, reorienting the orbital angular momentum vector, altering the orbital eccentricity, and even trading companions. However, from Figure 6, note that most of the mass accreted by a BH occurs when the orbital separation shrinks to ≲10 au, with a much longer interaction time.
• 3.  
  The accretion flow studied here is initially laminar, while AGN disk models typically assumed some α viscosity that parameterizes the role of isotropic, turbulent eddies that permeate the AGN disk. The characteristic length-scale of an eddy is ∼α H, and depending on the value of ${{ \mathcal M }}_{{\rm{w}}}$, this can be comparable to R_b. If so, the evolution of the embedded binary could be driven more strongly by the accretion of eddies with randomly distributed angular momenta than by the angular momentum supplied from the AGN disk velocity profile. If this is the case, mass accretion might occur similarly, but the spin evolution of the binaries would be stochastic. This is, however, a different flow geometry than that studied here and deserves attention in its own right.

While the relative importance of AGN disks as a formation scenario for BBH mergers may remain unresolved even after future LIGO/Virgo observing runs, the space-based gravitational-wave detector LISA may be sensitive to BBHs in AGN, a prospect that we will consider now. Although many different factors can impact the overall rate of systems observed by LISA, at sufficiently small orbital separations where gravitational-wave radiation dominates the orbital evolution of point masses such as BBHs, general relativity makes a clear prediction for the distribution of BBH orbital separations, P(a); since the strength of gravitational-wave radiation increases with decreasing orbital separation, ${(\dot{a}/a)}_{\mathrm{GR}}\sim {a}^{-4}$ (Peters 1964), a population of equal-mass BBH binaries will have P(a) ∼ a³. However, at sufficiently large separations, the gas drag forces will dominate. Taking p ≈ − 1, Equation (13) indicates that ${(\dot{a}/a)}_{\mathrm{drag}}\propto {a}^{5/2}$, when a ≪ R_b, producing an orbital separation distribution P(a) ∼ a^−7/2. By equating these two forces, we can find the critical orbital separation, a_crit, indicating the transition between the two regimes:

$\begin{eqnarray}\begin{array}{rcl}{a}_{\mathrm{crit}} & \simeq & \,25.7{\left(\displaystyle \frac{{m}_{1}}{30\,{M}_{\odot }}\right)}^{9/13}{\left({c}_{\infty }/20\,\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{-4/13}\\ & & \times {({\rho }_{\infty }/{10}^{-11}\,{\rm{g}}\,{\mathrm{cm}}^{-3})}^{-2/13}.\end{array}\end{eqnarray} \tag{ 18 }$

This corresponds to an orbital period ≈days or a gravitational-wave frequency ≈μHz. This is, unfortunately, well below the LISA sensitivity curve, which extends down to ≈mHz. Additionally, the strong dependence on a of both the gravitational-wave inspiral rate and our drag rate makes the dependence of a_crit on the gas density and sound speed extremely weak. The binaries may also grow significantly by accretion during their inspiral, which would cause a_crit to increase further. Equation (18) also assumes a purely hydrodynamic drag rate, while in reality the accretion flow is highly super-Eddington, which likely further decreases the efficiency of drag. Thus, if our a scalings predicted in Equation (13) hold, we find it unlikely for LISA to be able to detect embedded BBHs.

In recent years, the possibility that BBH mergers may be accompanied by EM signatures has intrigued the astrophysics community. There have been many proposed mechanisms for producing a signature, including the rapid accretion of relic disks post-merger (Perna et al. 2016; Schrøder et al. 2018; S. Schrøder et al. 2021, in preparation), the shock-heating of circumbinary disks due to the post-merger recoil (Corrales et al. 2010; Mink & King 2017), and mergers following a single-star progenitor (D'Orazio & Loeb 2018). In Bartos et al. (2017), they also studied the fates of embedded BBHs, and argued that these binaries could produce observable signatures by reaching super-Eddington luminosities via relativistic, beamed outflows. Interest in this scenario has been reinvigorated by the claimed association of an AGN flare with BBH merger GW190521. While this claim is tenuous (Ashton et al. 2020), it highlights one of the difficulties of this scenario: even if an embedded BBH produces a luminous signal, it must be distinguished from the dominant AGN emission. We note the recent work, Graham et al. (2020), which compared AGN flaring activity with BBH merger candidates observed by LIGO/Virgo in O3, and found 9 possible associations. These associations depend on the nature of EM signatures from merging embedded BBHs, the definition of what typical AGN flaring is, and the large localization windows from gravitational-wave observations. Given these considerations, we expect that the identification of EM signatures from embedded BBH mergers will only be made robust when we better understand what an associated flare would look like.

Radiation may be produced by an embedded BBH either by accretion during the inspiral phase or transiently immediately following the merger. We turn our attention, first, to the case of a steady-state luminosity. In many cases, the embedded BBH is supplied ⁷ gas at rates near the Eddington limit or significantly above it (as pointed out by other authors, e.g., Bartos et al. 2017; Stone et al. 2017). However, even in cases where the mass accretion rate can exceed the Eddington limit, the luminosity itself will still be limited because radiation becomes trapped within the flow. In this case, the accretion flow will be a geometrically and optically thick advection-dominated accretion flow ("ADAF"; Narayan & Yi 1994), producing an unusually soft blackbody spectrum. For an Eddington-limited pair of stellar-mass BHs to be observable, they must compete with the host AGN. The relative luminosity of the BBH accretor to the AGN is

$\begin{eqnarray}&&\displaystyle \frac{{L}_{\mathrm{BBH}}}{{L}_{\mathrm{AGN}}}\sim {10}^{-7}\left(\displaystyle \frac{{\lambda }_{\mathrm{Edd}}^{(\mathrm{SMBH})}}{0.1}\right)\left(\displaystyle \frac{M}{60\,{M}_{\odot }}\right){\left(\displaystyle \frac{{M}_{\mathrm{SMBH}}}{{10}^{8}\,{M}_{\odot }}\right)}^{-1}.\end{eqnarray} \tag{ 19 }$

For typical values of ${\lambda }_{\mathrm{Edd}}^{(\mathrm{SMBH})}$, M_BBH, and M_SMBH, this value is extremely small, and the accreting BBH will be indistinguishable from the host AGN unless the two are resolved.

Clearly, if the EM signatures from an embedded BBH are to be observed, it must shine brighter than the Eddington luminosity. The most natural way for this to happen is by producing a collimated, relativistic jet that has its emission beamed toward us (e.g., MacLeod et al. 2014). We can estimate an upper limit on the jet luminosity by assuming the flow is magnetically arrested, producing powerful Blandford–Znajek jets due to large poloidal magnetic fields (Blandford & Znajek 1977; Tchekhovskoy et al. 2011). Assuming this, we take

$\begin{eqnarray}&&{L}_{\mathrm{jet}}={\eta }_{\mathrm{MAD}}\dot{M}{c}^{2}\end{eqnarray} \tag{ 20 }$

where η_MAD is the jet efficiency for a magnetically arrested disk (MAD). The jet efficiency in MADs commonly reaches ∼200%–300% by tapping the rotational energy of the BH (Tchekhovskoy et al. 2011), so we assume a characteristic η_MAD = 2. Assuming that the radiative efficiency of the accretion flow is η ≈ 0.1, then the total jet luminosity is

$\begin{eqnarray}&&{L}_{\mathrm{jet}}\approx 20\left(\displaystyle \frac{{\eta }_{\mathrm{MAD}}}{2}\right){\left(\displaystyle \frac{\eta }{0.1}\right)}^{-1}{L}_{\mathrm{Edd}}.\end{eqnarray} \tag{ 21 }$

If this jet is directed toward us, the emission will be beamed within an opening angle θ ∼ 1/Γ, where Γ is the Lorentz factor of the jet. If we measure the flux from this jet at Earth, then the inferred isotropic luminosity is

$\begin{eqnarray}&&{L}_{\mathrm{iso}}\sim \displaystyle \frac{4\pi }{{\theta }^{2}}\displaystyle \frac{{L}_{\mathrm{jet}}}{2}={{\rm{\Gamma }}}^{2}2\pi {L}_{\mathrm{jet}}.\end{eqnarray} \tag{ 22 }$

If we compare this boosted luminosity to the host AGN luminosity, then we find that

$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{L}_{\mathrm{iso}}}{{L}_{\mathrm{AGN}}} & \sim & {10}^{-5}{{\rm{\Gamma }}}^{2}\left(\displaystyle \frac{2}{{\eta }_{\mathrm{MAD}}}\right){\left(\displaystyle \frac{0.1}{\eta }\right)}^{-1}\left(\displaystyle \frac{{\lambda }_{\mathrm{Edd}}^{(\mathrm{SMBH})}}{0.1}\right)\\ & & \times \left(\displaystyle \frac{M}{60\,{M}_{\odot }}\right){\left(\displaystyle \frac{{M}_{\mathrm{SMBH}}}{{10}^{8}\,{M}_{\odot }}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 23 }$

If we require that L_iso ≳ L_AGN, then the corresponding Γ for which this is achieved is Γ ≈ 316. This is in the upper range of Γ values associated with gamma-ray bursts (Lithwick & Sari 2001; Gehrels et al. 2009), and corresponds to an opening angle θ ∼ 018. Even if embedded BBHs could produce jets with this Lorentz factor, we would only see ×θ²/4 ∼ 0.0002% of them due to the inclination-dependence of the emission. Furthermore, while Bartos et al. (2017) suggest that embedded BBHs form gaps, allowing jets to escape unimpeded, Figure 2 suggests most binaries will not form gaps, and Figure 6 suggests most binaries merge in the outer regions (≳0.5–1 pc) of the disk. In these regions, ${{ \mathcal M }}_{{\rm{w}}}\sim 0.1$, and H is roughly 2 orders of magnitude larger than the binary sphere of influence. This provides a large column density of material for the jet to traverse, likely resulting in significant mass-loading, which would lower the resulting Γ. This suggests that it is difficult to detect a jet associated with an embedded BBH unless it occurs at frequencies separate from the dominant AGN emission. We note that the EM signatures associated with these jets have been studied in detail by two recent papers: Zhu et al. (2021), for the case of embedded neutron stars, and Perna et al. (2021), for jet propagation from both BHs and neutron stars. In both works, they emphasize the issue of distinguishing the transient emission from the AGN emission, and find that the transient can more easily outshine the AGN in the X-ray, infrared, and radio bands.

Finally, we consider the possibility of transient EM signatures that occur during or immediately following the merger. There are two possibilities; the first possibility is that the binary experiences a transient spike in its accretion rate and thus luminosity post-merger (Milosavljević & Phinney 2005). However, these studies were focused on thin circumbinary disks with carved out cavities, and in the thick flows that accompany these binaries, no such spike is likely to exist. We therefore ignore this possibility. The other possibility is that the sudden decrease in gravitational potential at the onset of the merger rapidly shocks the surrounding gaseous envelope, producing an intense photospheric emission.

The sudden mass loss accompanying the BBH merger is associated with a drop of potential energy in the surrounding high-density envelope,

$\begin{eqnarray}&&{\rm{\Delta }}E={\int }_{{R}_{\mathrm{bw}}}\displaystyle \frac{G{\rm{\Delta }}{M}_{\mathrm{merger}}}{r}\rho (r){dV}.\end{eqnarray} \tag{ 24 }$

We assume that within the BBH sphere of influence, characterized by R_bw (Equation (11)), the density profile is Bondi-like, i.e., $\rho {(r)\sim {\rho }_{\infty }(r/{R}_{\mathrm{bw}})}^{-3/2}$ (Shapiro & Teukolsky 1983). This yields a change in potential energy,

$\begin{eqnarray}&&{\rm{\Delta }}E=4\pi {\rho }_{\infty }G{\rm{\Delta }}{M}_{\mathrm{merger}}{R}_{\mathrm{bw}}^{2}.\end{eqnarray} \tag{ 25 }$

If we assume the flow is Bondi-like in the region of interest and that roughly ∼50% of the BBH mass is lost upon merger, then this change in energy can be rewritten as

$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}E & \approx & 8.1\times {10}^{40}\,\mathrm{erg}{\left(\displaystyle \frac{M}{60\,{M}_{\odot }}\right)}^{3}\left(\displaystyle \frac{{\rho }_{\infty }}{{10}^{-18}\,{\rm{g}}\,{\mathrm{cm}}^{-3}}\right)\\ & & \times {\left(\displaystyle \frac{{c}_{\infty }}{25\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)}^{-4}.\end{array}\end{eqnarray} \tag{ 26 }$

If we crudely assume that this energy is lost within a sound-crossing time Δt = R_b/c_∞, then the luminosity produced by the shock (≈ΔE/Δt) is

$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{shock}} & \approx & 1.6\times {10}^{32}\,\mathrm{erg}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{M}{60\,{M}_{\odot }}\right)}^{2}\left(\displaystyle \frac{{\rho }_{\infty }}{{10}^{-18}\,{\rm{g}}\,{\mathrm{cm}}^{-3}}\right)\\ & & \times {\left(\displaystyle \frac{{c}_{\infty }}{25\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 27 }$

Since a near-Eddington, 10⁸ M_⊙ host SMBH should have a luminosity of roughly L_AGN ≈ 10⁴⁵–10⁴⁶ erg s⁻¹, it is implausible for the luminosity of the post-merger shock to exceed that of the AGN.

One of the uncertainties in the accretion flow embedding the binaries is the role of feedback. Feedback can be deposited into the gas supply either radiatively or mechanically, through winds and relativistic jets, and can effect the mass accretion rate and efficiency of drag.

If an accreting binary produces a jet, it is possible for the jet to carve out a cavity and suppress the gas supply (e.g., Ramirez-Ruiz et al. 2002). If the force exerted on the ambient medium by a relativistic jet ∼L_jet/c, then it will exert a ram pressure ∼L_jet/c θ² r², where θ is the opening angle of the jet. By assuming that, as the jet plunges through the medium, it becomes subrelativistic and spreads laterally, we can take θ² = 4π. We can set the jet ram pressure equal to the gas pressure of the ambient medium to determine the radius of the jet-medium interface (denoted r_j ),

$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{L}_{\mathrm{jet}}}{4\pi {r}_{j}^{2}c} & = & {\rho }_{\infty }{c}_{\infty }^{2},\\ {r}_{j} & = & \sqrt{\displaystyle \frac{{L}_{\mathrm{jet}}}{4\pi {\rho }_{\infty }{c}_{\infty }^{2}c}}.\end{array}\end{eqnarray} \tag{ 28 }$

To compare this to the relevant physical scale, we can write r_j in terms of the Bondi radius and physical values,

$\begin{eqnarray}\begin{array}{rcl}{r}_{j}/{R}_{{\rm{b}}} & \approx & 56{\left(\displaystyle \frac{{\lambda }_{\mathrm{Edd}}^{(\mathrm{jet})}}{20}\right)}^{1/2}{\left(\displaystyle \frac{M}{60\,{M}_{\odot }}\right)}^{-1/2}\\ & & \times {\left(\displaystyle \frac{{\rho }_{\infty }}{{10}^{-18}\,{\rm{g}}\,{\mathrm{cm}}^{-3}}\right)}^{-1/2}\left(\displaystyle \frac{{c}_{\infty }}{25\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right).\end{array}\end{eqnarray} \tag{ 29 }$

In this relation, we have expressed the jet luminosity in terms of its Eddington ratio, ${\lambda }_{\mathrm{Edd}}^{(\mathrm{jet})}$. For the assumed ambient conditions (typical of the AGN disk at ≈1 pc) and Eddington ratio (for an upper limit, see Equation (21)), r_j is much larger than R_b. This suggests that if a powerful jet is produced by the binary, it should easily carve out a large cavity and impede the gas supply. This could at least transiently diminish the accretion rate and may impose a duty cycle on the growth of the BBH. Alternatively, we can use r_j < R_B as the condition required for a cavity not to form. Then, the maximum Eddington ratio of the jet is

$\begin{eqnarray}&&{\lambda }_{\mathrm{Edd}}^{(\mathrm{jet})}\lt 0.0064\left(\displaystyle \frac{M}{60\,{M}_{\odot }}\right)\left(\displaystyle \frac{{\rho }_{\infty }}{{10}^{-18}\,{\rm{g}}\,{\mathrm{cm}}^{-3}}\right){\left(\displaystyle \frac{{c}_{\infty }}{25\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)}^{-2}.\end{eqnarray} \tag{ 30 }$

This is an extremely low jet efficiency and would likely require a BH with a very small spin (Narayan et al. 2003). However, closer to ≈10⁻¹–10⁻² pc, the gas density can reach ≈10⁻¹⁴ g cm⁻³ (Sirko & Goodman 2003), which would result in a minimum Eddington ratio of 64. This suggests that the ability of a jet to impede the gas supply strongly depends on where the BBH is in the disk; at larger distances, it likely induces a duty cycle in the accretion flow, while at smaller distances the jet will be quickly extinguished.

The results of this work, particularly the evolutionary tracks laid out in Section 4.2, have uncertainties that deserve emphasis. Here are a few essential ingredients required to properly understand the evolution of embedded BBHs:

• 1.  
  Understanding how drag and mass accretion rates change in the super-Eddington regime is critical to understanding the evolution of embedded BBHs. The ratio of the accretion to inspiral timescale is what determines the final mass of the binary upon merger and will be altered if radiation is realistically taken into account. Additionally, the radiative efficiency affects both the mass accretion rate and the evolutionary timescale. While we use a constant radiative efficiency of η = 0.01, it will be higher in the sub-Eddington regime, but can be even lower in the highly super-Eddington regime. While the details of the evolutionary tracks depend on how the drag in the super-Eddington regime is prescribed, two main features are clear: that embedded BBHs can grow significantly, possibly entering the PISN mass gap before merging, and that the time to merger is ≲1–30 Myr.
• 2.  
  The true nature of AGN disks, particularly in the outer Toomre-unstable region, remains an open question. While we favor the profiles of Sirko & Goodman (2003), they use an unspecified pressure source to stabilize this region. This provides a useful model for making predictions but remains a crude estimate, which should be explored further. Additionally, as discussed in Section 4.2, the time-dependence of AGN disks is a major uncertainty. As found in the simulations of Angles-Alcazar et al. (2020), at parsec scales, the disk can reorient itself on ≈Myr timescales, which may alter the occupation fraction of BHs in the AGN disk or limit the total amount of time the binaries have to evolve in the disk.
• 3.  
  We assumed that the embedded BBH is allowed to evolve through the AGN undisturbed, but in reality it dynamically couples to the surrounding stellar population. Three-body scattering may merge the BBH faster (Stone et al. 2017) or may reorient the binaries with respect to the surrounding flow, which can alter the mass and spin evolution of the binaries prior to merger. Additionally, as discussed in Section 4.1, a proper understanding of the BH distribution near the SMBH will allow us to better understand how many BHs can be captured dynamically by the AGN and where their initial position in the disk is, which can strongly affect their subsequent evolution.
• 4.  
  We have neglected multiple effects in our hydrodynamic model. Importantly, these include the noninertial forces of coriolis and centrifugal forces, which will likely modify the flow. Additionally, properly modeling the vertical stratification of the flow is also likely to be important. Even given a more realistic domain, there are large uncertainties pertaining both to the size of the sink radius (which will affect the convergence of the mass accretion rate, e.g., Appendix A) and on the neglected physics. In particular, a proper treatment of the thermodynamics (i.e., via a cooling prescription or better yet radiation transport) will affect the mass accretion and drag rates (e.g., Li et al. 2022), and the feedback processes (such as relativistic jets or outflows launched on subgrid scales) may also alter the flow substantially.

To summarize, we describe the chronology of embedded BBHs as they evolve hydrodynamically through the host AGN disk:

• 1.  
  Birth. A pair of BHs either are born in the disk (typically ≳0.1–1 pc) or are dynamically captured. If the NC harbors a BH cusp, it is possible for the BHs to be initially captured at much smaller radii. As discussed in Section 4.1, it is unclear whether or not a BH cusp will have been formed during the AGN phase of the SMBH, particularly for higher-mass SMBHs where the relaxation timescale of the surrounding stellar population can be longer than the age of the universe. For this reason, we find it more likely that BBHs first become embedded in the outer regions of the AGN disk.
• 2.  
  Accretion. As embedded BBH orbits prograde with the AGN disk, it will become engulfed in a gaseous wind that has an asymmetric velocity profile in the rest frame of the BBH. This is depicted in Figures 1, 2, and 3, where the characteristics of the accretion flow are primarily defined by ${{ \mathcal M }}_{{\rm{w}}}$, the Mach number of the wind in the BBH rest frame. At large ${{ \mathcal M }}_{{\rm{w}}}$, the ram pressure of the wind dominates the accretion flow, and at small ${{ \mathcal M }}_{{\rm{w}}}$, the accretion flow becomes Bondi-like. In the outer, Toomre unstable regions of an AGN disk, ${{ \mathcal M }}_{{\rm{w}}}$ is often small enough that Bondi accretion provides an adequate approximation (Figure 2). When ${{ \mathcal M }}_{{\rm{w}}}$ is large enough (≫1), the accretion flow becomes planar, likely resembling a pair of thin mini-disks and potentially opening an annular gap in the AGN disk (Baruteau et al. 2011; Li et al. 2021). Our analysis is invalid in this regime, which likely holds for AGN disks that are geometrically thin and are hosted by lighter (∼10⁶ M_⊙) SMBHs. This makes our results most salient for binaries embedded within the disks of heavier (≳10⁷ M_⊙) SMBHs. For most regions in the disk, a BBH will accrete at super-Eddington rates (Section 4.2); if the BBH accretes at the Eddington limit for ≳5 Myr, it can double its mass.
• 3.  
  Inspiral. The inspiral of the BBH is most strongly governed by mass accretion and drag, with the drag generally contributing more than the mass accretion (Figure 4). In Section 3.3, we provide analytic expressions with best-fit numerical coefficients that reasonably reproduce the inspiral rate measured in our hydrodynamic simulations. We find that typical BBHs, beginning their journey anywhere between 0.1 and 2 pc in the AGN disk, will merge within ≲1–30 Myr.
• 4.  
  Migration. Concurrent to accretion and inspiral, the embedded BBHs will migrate in the AGN disk (typically inwards). For the initial distances we considered (0.5–2 pc), we found that our BBHs migrated only marginally before merging.
• 5.  
  Multimessenger signatures. In Section 4, we found that BBHs can often double their mass during their inspiral. This can alter the resulting BH mass spectrum detected by LIGO, allowing embedded BBHs to enter the PISN mass gap, as discussed in Section 5.1. In Section 5.2, we studied the possible EM signatures accompanying an accreting BBH or occurring transiently post-merger. In general, it is difficult for the BBH to outshine its host AGN. If the BBH produced a relativistic jet, it would have to have a large Lorentz factor to be detectable, leaving only a very small fraction of these sources detectable due to the inclination-dependence of the beamed emission. Shock-heating in the AGN disk due to the post-merger recoil can produce EM transients, but with luminosities far less than that of the AGN.