Orbital evolution of binary black holes in active galactic nucleus disks: a disk channel for binary black hole mergers?
DDraft version January 26, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Orbital evolution of binary black holes in active galactic nucleus disks: a disk channel for binary blackhole mergers?
Ya-Ping Li, Adam M. Dempsey, Shengtai Li, Hui Li, and Jiaru Li Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA (Dated: January 26, 2021, submitted to ApJ)
ABSTRACTWe perform a series of high-resolution 2D hydrodynamical simulations of equal-mass binary blackholes (BBHs) embedded in active galactic nucleus (AGN) accretion disks to study whether thesebinaries can be driven to merger by the surrounding gas. We find that the gravitational softeningadopted for the BBH has a profound impact on this result. When the softening is less than tenpercent of the binary separation, we show that, in agreement with recent simulations of isolated equal-mass binaries, prograde BBHs expand in time rather than contract. Eventually, however, the binaryseparation becomes large enough that the tidal force of the central AGN disrupts them. Only whenthe softening is relatively large do we find that prograde BBHs harden. We determine through detailedanalysis of the binary torque, that this dichotomy is due to a loss of spiral structure in the circum-singledisks orbiting each BH when the softening is a significant fraction of the binary separation. Properlyresolving these spirals – both with high resolution and small softening – results in a significant sourceof binary angular momentum. Only for retrograde BBHs do we find consistent hardening, regardlessof softening, as these BBHs lack the important spiral structure in their circum-single disks. Thissuggests that the gas-driven inspiral of retrograde binaries can produce a population of compact BBHsin the gravitational-wave-emitting regime in AGN disks, which may contribute a large fraction to theobserved BBH mergers.
Keywords: accretion disks—galaxies: active – gravitational waves – stars: black holes—planet-diskinteractions INTRODUCTIONSince the first discovery of a transient gravitational-wave signal by the two detectors of the Laser Interfer-ometer Gravitational-Wave Observatory (LIGO;Abbottet al. 2016), more than fifty mergers have been dis-covered during LIGO/Virgo’s first three observing runs(Abbott et al. 2019, 2020a). Active galactic nucleus(AGN) disks have been proposed as promising loca-tions for producing some of the detected stellar massbinary black hole (BBH) mergers (McKernan et al. 2012,2014, 2018, 2019; Bartos et al. 2017b; Stone et al. 2017;Leigh et al. 2018; Secunda et al. 2019, 2020; Yang et al.2019a,b; Gr¨obner et al. 2020; Ishibashi & Gr¨obner 2020;Tagawa et al. 2020a,b; Abbott et al. 2020b). Other fa-
Corresponding author: Ya-Ping [email protected] vored scenarios include isolated binary evolution in abinary star system (e.g., Belczynski et al. 2010; Man-del & de Mink 2016; Belczynski et al. 2016), dynamicalevolution of triple or quadruple systems (e.g., Antoniniet al. 2017; Liu & Lai 2017, 2019; Michaely & Perets2019; Fragione & Kocsis 2019), and dynamical forma-tion in which the black holes (BHs) undergo a chanceencounter in a dense stellar environment such as galacticnuclei, globular clusters or open clusters (e.g., O’Learyet al. 2009; Wang et al. 2016; Banerjee 2017; Rodriguezet al. 2018).An AGN disk is a favorable location for BBH mergersbecause the surrounding gaseous disks may harden ex-isting BBHs (McKernan et al. 2014; Bartos et al. 2017b;Stone et al. 2017; Yang et al. 2019a). AGN disk-assistedBH mergers have distinct properties that could differen-tiate them from other merger channels. These includeheavier BBH mergers (Yang et al. 2019a; Abbott et al.2020b), a spatial correlation with AGN host galaxies a r X i v : . [ a s t r o - ph . H E ] J a n Li et al. (Bartos et al. 2017a), large spin magnitudes (McKernanet al. 2012), and possible electromagnetic counterpartsdue to the BHs accreting from the surrounding gaseousdisk (Stone et al. 2017; Bartos et al. 2017b; McKernanet al. 2019; Graham et al. 2020).The formation of BBHs within the AGN disk can bedue to migration traps (McKernan et al. 2012). As BHsaccumulate near these traps, they can form BBH sys-tems (Bellovary et al. 2016; Secunda et al. 2019). Inaddition, BBHs can be formed via isolated evolution ofthe stellar binary systems (Mandel & de Mink 2016), orformed dynamically through the chance encounters ofBHs in galactic nuclei (O’Leary et al. 2009), or formedwithin the self-gravitating disk itself (Stone et al. 2017).Most works on the AGN-assisted channel adopt semi-analytic results for the drag force in a uniform gaseousbackground medium (Kim & Kim 2007; Kim et al.2008) or from a circum-binary disk (CBD) (Goldreich &Tremaine 1980; Cuadra et al. 2009; Roedig et al. 2012;D’Orazio et al. 2013) to conclude that gas surroundingthe BBH facilitates its merger. This requires confirma-tion from high resolution numerical simulations since thereal situation for an embedded binary in an AGN diskis much more complicated.The expectation of BBH coalescence stems from earlyhydrodynamical simulations of isolated binaries (e.g.,MacFadyen & Milosavljevi´c 2008; Shi et al. 2012). How-ever, these simulations did not include the binary it-self in the computational domain and so missed the im-portant source of angular momentum coming from thecircum-single disks (CSDs). Later work showed that in-cluding and appropriately resolving CSDs in the simula-tion causes isolated binaries to expand rather than con-tract (Tang et al. 2017; Mu˜noz et al. 2019, 2020; Duffellet al. 2020). Similar results have been obtained from 3Dhydrodynamical simulations (Moody et al. 2019). More-over, early simulations of isolated binaries were typicallyevolved for a short amount of time. This was shown tocause transient hardening whereas in steady-state thebinary expands (Miranda et al. 2017).Recently, there have been some isolated binary simu-lations which showed that the binary can still be hard-ened by the gaseous disk with a small binary mass ratio(Duffell et al. 2020; Derdzinski et al. 2021) or a cold disk(i.e., small disk scale height, Tiede et al. 2020; Heath& Nixon 2020). Chen et al. (2020a) found analyticallythat the radiation field from AGN also causes the or-bits of the accreting super-massive black hole binaries(SMBHBs) with a small mass ratio to shrink by the“Poynting-Robertson” drag effect.To our knowledge, there has only been one work thathas used 2D hydrodynamical simulations to study an embedded binary star in an AGN disk (Baruteau et al.2011). While the outcome of their simulations agreedwith the semi-analytic and early isolated binary results,the adopted resolution and gravitational softening couldnot adequately resolve the innermost regions (i.e., CSDregions) around each star. Therefore, the AGN diskchannel for BBH and SMBHB mergers remains an activeresearch topic.In this work, we have performed a series of high-resolution 2D hydrodynamical simulations of BBHs em-bedded in thin accretion disks. The resolution and gravi-tational softening in these simulations are chosen so thatthe CSDs around each BH are properly resolved. Ourprimary goal is to examine whether or under what con-ditions a BBH will harden in an AGN disk. Motivatedby the isolated binary simulations, we quantitatively ex-plore the effects of the CSD structures on the binaryevolution.The paper is organized as follows. In Section 2 we in-troduce the numerical setup of our simulations, and inSection 3 we present our numerical results. We explorethe binary evolution for different gravitational soften-ing scales and prograde/retrograde orbits. Finally, wesummarize the main finding of our paper, and then dis-cuss the limitations of our present work and the futuredirection in Section 4. METHODWe numerically solve the 2D continuity and Navier-Stokes equations using the hydrodynamical code
LA-COMPASS (Li et al. 2005, 2009), to study a BBH em-bedded in a gaseous disk that surrounds a super-massiveblack hole (SMBH). For simplicity, only one BBH is in-cluded in the disk.In our simulations, we use a geometrically thin andnon-self-gravitating disk. We choose a 2D cylindricalcoordinate system ( r, ϕ ), with the origin located at theposition of the central SMBH.2.1.
Initial Conditions
We adopt a simple model in which the disk tempera-ture follows T disk ∝ r − . With this, the disk aspect ratio h ( ≡ H/r ) is constant throughout the whole disk andwe choose h = 0 .
08. Here H is the disk scale height.The disk surface density has an initial profile ofΣ( r ) = Σ (cid:18) rr (cid:19) − . , (1)where Σ is the gas surface density at r , and r is theinitial radius of the binary circular orbit for the BBHcenter of mass (COM). The 2D velocity vector of thegas is v = ( v r , v ϕ ), and angular velocity is Ω = v ϕ /r . BH in AGN Disks α -prescription for the gas kinematic vis-cosity ν = αH Ω (Shakura & Sunyaev 1973). Such acombination of Σ, T disk , and ν ensures an initial steadystate of disk accretion.We choose the unit of length to be r , and the unitof mass to be the SMBH mass, M smbh . We set Σ =10 − M smbh /r , which corresponds to a total disk massof 10 − M smbh . With these values, the gas self-gravityis negligible. 2.2. Gravitational Interaction
We initialize an equal mass BBH on a circular or-bit with a semi-major axis of a = κR H , where R H =(2 q/ / r is the BBH Hill radius, and q is the massratio of one BH to the SMBH. We keep κ <
1, so thatthe SMBH does not strongly affect the gas dynamicsin the vicinity of the BBH. The binary orbital period is P bin = (cid:112) κ / P com in units of the binary’s COM orbitalperiod around the SMBH, P com .To the gas momentum equation, we add the gravita-tional potential, φ = − GM smbh | r | + (cid:88) i =1 (cid:34) − GqM smbh ( | r bh ,i − r | + (cid:15) ) / + q Ω ,i r bh ,i · r (cid:35) . (2)The location and orbital frequency of each BH with re-spect to the SMBH are given by r bh ,i and Ω bh ,i , where i = 1 , G is thegravitational constant. The first term is the potentialdue to the central SMBH, while the first term in thebracket is the direct potential from each BH. The sec-ond term in the bracket is the indirect potential aris-ing from our choice of coordinate system. We apply agravitational softening, with length scale (cid:15) , to each BHpotential. Based on the isolated binary simulations, asmall gravitational softening scale (i.e., (cid:28) a ) is neededto resolve the CSD in the vicinity of each BH. Here, wemainly focus on a softening of 0 . a , which is compara-ble to the softening scale in isolated binary simulations(e.g., Mu˜noz et al. 2019, 2020; Tiede et al. 2020). Tobetter understand the influence of the softening on ourresults, we also explore a range of values.Applying such a small softening requires a super reso-lution around each BH to resolve their CSDs. For typicalstellar mass BH and SMBH masses, q (cid:39) − ∼ − and R H (cid:39) . ∼ . r . Given that the size ofthe CSD region is much smaller than R H , resolving thisregion for realistic parameters is unfeasible. Becauseof this, we adopt a larger mass ratio of q = 2 × − ,disk viscosity α = 6 . × − , and disk scale height h = 0 .
08. This setup maintains the same gap pro-file with the more realistic system of q = 4 × − , α = 8 × − and h = 0 .
01 (Lin & Papaloizou 1993;Crida et al. 2006; Baruteau et al. 2011).2.3.
Gas Accretion
Our BH accretion prescription gradually removesmass and momentum from the gas within a distanceof r acc = 0 . r to each BH, which is comparable tothe BH softening scale. Material that is accreted by aBH is not added to its mass or momentum. We fix thetime scale of this removal to τ − ,i = 5 Ω bh ,i unless oth-erwise stated. We have tested that this choice of τ acc ,i – even τ acc ,i → ∞ – does not change our main resultssignificantly (see also Li et al. 2021 for the effect of τ acc on the accretion rate). We discuss other removal ratesin Section 3.5.2.4. Resolution and Boundary Conditions
We resolve the disk with a uniform radial grid of n r = 2048 points in a radial domain between [0 . , . r ,and a uniform azimuthal grid with n φ = 8192 points.With such a high resolution, the BBH Hill radius canbe resolved with ∼
140 cells. We also test that a largerradial domain with r ∼ [0 . , . r in a logarithmic ra-dial grid of n r = 3072 does not have a significant impacton the binary dynamics.At both boundaries we fix the gas to its initial ax-isymmetric configuration. To facilitate this, we gradu-ally remove any waves near both boundaries with wavekilling zones (de Val-Borro et al. 2006).Our simulations are performed in two steps. The BBHis first held on a fixed circular orbit for 500 orbits, sothat the disk can reach a quasi-steady state. Duringthis time, the BBH opens a gap of similar depth andshape to an equivalent mass single planet. After 500orbits, the force from the disk on the BBH is includedand the binary can dynamically evolve. For the entiretyof the simulation, we do not allow the time step to ex-ceed 1 /
400 of the binary’s current internal orbital pe-riod. This ensures that the dynamics in the CSDs isadequately resolved in time (Cresswell & Nelson 2006;Baruteau et al. 2011). NUMERICAL RESULTSTable 1 provides an overview of all of our simulations.We first study the cases where the binary and disk an-gular momenta are aligned. This is the prograde caseshown in Table 1.3.1.
Gas Structures
As a typical case, we choose Model A to show the gassurface density in Figure 1. On scales larger than R H , Li et al.
Figure 1.
Snapshots for the gas surface density in Model A at 500 COM orbits. From left to right, each panel zooms incloser to the BBH. The dashed white circle shows the BBH Hill radius. Overlaid arrows show the gas velocity streamline in theco-moving frame of the BBH COM. The dashed red circles in the right panel correspond to the softening scale for each BH.
Table 1. Model parameters and binary evolution
Model a (cid:15) accretion prograde/ Comments( r ) ( a ) retrogradeA 0.02 0.08 on prograde outspiralingAr 0.02 0.08 on a prograde outspiralingB 0.02 0.4 on prograde inspiralingC 0.02 0.2 on prograde outspiralingD 0.02 0.08 off prograde outspiralingE 0.04 0.04 b on prograde outspiralingF 0.04 0.2 b on prograde inspiralingG 0.04 0.04 b on retrograde inspiralingH 0.04 0.2 b on retrograde inspiraling Note — a A smaller removal rate of 0 . bh ,i and a largerradial domain of r ∼ [0 . ,
4] for the BBH is adopted to testthe effect of binary accretion. b We keep the same physical length of gravitational softeningas Model A or B. For the first five Models (A, Ar, B, C & D), a = 0 . R H , while for the remaining Models a = 0 . R H . the disk is nearly indistinguishable from a single largemass planet (cf. Dempsey et al. 2020). On smaller scales,the gas streamlines reveal that gas circulates around theBBH COM. Based on this, there is a clear “CBD” regionjust inside R H . This is different from the viscously con-trolled CBD in isolated binary simulations (e.g., Tanget al. 2017; Mu˜noz et al. 2019), as the existence of thecentral SMBH drives large-scale spiral arms that feedthe CBD.Zooming even closer in to the binary shows clearCSDs. With such a small softening for Model A, theCSD around each BH is properly resolved by our highresolution, and contains several spiral arms. In particu-lar, the inter-spiral arm connecting the two BHs is quiteprominent, as are the two trailing spiral arms that con- nect to the large scale spirals through the CBD region.There is not a sharp cavity around the CSD region asshown in many isolated binary simulations (e.g., Tanget al. 2017; Mu˜noz et al. 2019), which is probably dueto the shallow depth of the gap carved by the BBH.3.2. Binary Dynamics
In Figure 2, we show the time evolution of the BH-BH and BBH-SMBH orbits for Model A. The orbit ofthe BH-BH binary is characterized by the semi-majoraxis a bin and eccentricity e bin , while the BBH-SMBHbinary has semi-major axis a com and is circular through-out its evolution. After release, a bin increases with timeas shown in the upper panel. This is consistent withmany recent isolated binary simulations (Miranda et al.2017; Moody et al. 2019; Mu˜noz et al. 2019, 2020; Duffellet al. 2020; Tiede et al. 2020). The shaded region showsthe instantaneous separation between the two BHs. Wecan see that after ∼ e bin ∼ . ∼ . a (Model B), the binary gradually inspirals atnear zero eccentricity until it reaches the softening scaleand stalls (as seen in e.g., Baruteau et al. 2011). Forthis case, the COM migration rate is similar to a single q = 4 × − planet (e.g., Chen et al. 2020b). Thissuggests that the global migration is mainly controlledby the torque far beyond the binary’s CBD region.To explore the effect of the binary release on the sepa-ration evolution, we show the long-term evolution rate of˙ a bin /a bin contributed by the gaseous disk in the fourthpanel of Figure 2, which is defined as˙ a bin a bin = − ˙ ε bin ε bin , (3) BH in AGN Disks Figure 2.
Binary evolution for Models A (black), and B(red). The first panel shows the time evolution of the BBHinstantaneous separation (shaded region) and semi-majoraxis (solid lines), while the second panel shows the BBHeccentricity. The third panel plots the time evolution of theBBH COM semi-major axis, and the last panel shows therate of change of the BBH semi-major axis due to the diskforces. The dashed line in the first panel corresponds to thesoftening scale in Model B. In Models A, the BBH expandswhile increasing its eccentricity, while in Model B the BBHcontracts to the softening scale. where ε bin = − GqM smbh /a bin is the specific energy ofthe binary, and dε bin dt = ( ˙ r bin , − ˙ r bh , ) · ( f grav , − f grav , ) . (4)Here, f grav , and f grav , are the gravitational force foreach BH from the disk. We have ignored the contribu-tion from accretion since we do not add the accretedmass and momentum to the binary. By the time thebinary is released, ˙ a bin is converged in time and there isno impulsive kick upon release.3.3. Gravitational Torque
In the following, we analyze the torque acting on thebinary to understand what causes the binary inspiral-
Figure 3.
Cumulative gravitational torque for the binaryevolution as a function of the distance to the binary’s COM δr for Models A, B, and C. The dashed lines correspond tothe three torquing regions discussed in the text: the inter-spiral arms, CSDs, and CBD. Gas outside of δr ∼ . a ,provides little to no torque on the BBH. The sign of T g for δr (cid:29) T g <
0) orexpands ( T g > ing/outspiraling for the different softening lengths. Todecompose the torque density spatially around the bi-nary, we follow the isolated binary literature and com-pute the gravitational torque density in three differentregions: inter-spiral arms ( δr ≤ . a ), CSDs (0 . a <δr ≤ a ), and CBD ( δr > a ). The torque densityfrom each computational cell acting on the binary is d T g /dA = ( r bh , − r bh , ) × ( f grav , − f grav , ).In Figure 3, we show the the cumulative radial torquewith a time averaging around 500 COM orbits, T g ( <δr ) = (cid:82) δr ( dT g /dA ) dA with δr = | r − r com | , for ModelsA,B, and C. All models, regardless of softening, find thatgas inside of 0 . a and outside of a torque the BBHdown, while gas between 0 . a and a torque the binaryup. However, the amount that each region contributesto the total torque is highly softening dependent. As thesoftening becomes smaller, the total amount of positivetorque coming from inside of a increases until it over-comes the amount of negative torque coming from theCBD region. The transition from CBD to CSD dom-inated angular momentum transfer, and consequentlyfrom inspiraling to outspiraling, occurs around a soften-ing of (cid:15) ≈ . a .To better understand the spatial contribution of thegravitational torque, we follow Mu˜noz et al. (2019) (seealso Tang et al. 2017; Tiede et al. 2020) to reconstruct Li et al.
Figure 4.
Stacked gas surface density (left panels) and specific gravitational torque density dT g /dA (middle and right panels)at t = 500 orbits for binary’s COM frame averaged over 100 snapshots. The coordinate ( x (cid:48) − x com , y (cid:48) − y com ) is the frameco-rotating with the BBH with its COM located at ( x com , y com ). The summation of dT g /dA (middle panels) over the verticalcoordinate gives the symmetric component of the torque density dT S g /dA (right panels), showing in the second and fourthquadrants. Upper panels: softening scale is 0 . a for each BH. Lower panels: softening scale is 0 . a . The white and blackdashed circles at δr = 0 . a , δr = a , and δr = 2 . a , in these plots delineate the different zones for the torque spatialdecomposition. The red dashed circles denote the corresponding gravitational softening scales. the 2D torque density distribution around 500 COM or-bits. We stack the gas surface density Σ and torquedensity dT g /dA in the binary internal co-rotating frameover one binary’s COM period with 100 snapshots, andshow the results in Figure 4. The upper and lower panelscorrespond to the small (Model A) and large (Model B)softening cases, respectively. The symmetric componentof the torque density with respect to the vertical coordi-nate, are shown in the right panels. The total torque onthe BBH can be obtained by summing all four quadrantsin the middle panels, or by summing the upper left andlower right quadrants in the right panels.From the stacked surface density shown in the leftpanels, the coherent inter-spiral arm and outer spiralarms connecting the two BHs are quite prominent. Incontrast, the spiral arms on the large scale (outside Hillradius of the binary indicated by outer white circles) aresmoothed out as they are not stationary in the binary’srotating frame. We find that the torque density can be well describedby four distinctive quadrants. It is positive in the firstand third quadrants, while negative in the second andfourth quadrants. By comparing with the stacked sur-face density plots, both the inter-spiral arms and theouter-spiral arms contribute negative to the total grav-itational torque, while the two-tips of the outer-spiralarms located in the second and fourth quadrants con-tribute positive torque.When the softening is small, the positive contribu-tion from the CSD region dominates over the negativecontribution from the spiral arms, as shown from thesymmetric torque density plot (the upper right panel ofFigure 4), and the radial cumulative torque shown inFigure 3. Inside of 0 . a , the interconnecting spiraltorques the binary down, while between 0 . a and a the CSD spiral arms add enough angular momentum tothe binary to overcome the negative torque from inside0 . a . Outside a , the CBD spirals do not contribute BH in AGN Disks
Binary Eccentricity and Destruction
An intriguing feature for the outspiraling binary is therapid increase in eccentricity after its semi-major axisexceeds ∼ . r . The maximum eccentricity duringthis stage can be as large as 0.5, as shown in the secondpanel of Figure 2. Since the semi-major axis of the bi-nary a bin does not increase significantly, the minimumseparation can be still be smaller than the initial valueof a .If we calculate ˙ e bin due to the disk, the rate we findis about one order of magnitude smaller than what isshown in Figure 2. This suggests that the rapid in-crease in binary eccentricity is not due to the disk,but to the interaction with the central SMBH. Indeed,there are BBH-destroying instabilities in hierarchical co-planar triple systems that occur at our chosen masseswhen a bin ∼ . a com ≈ . R H (Eggleton & Kisel-eva 1995; Mardling & Aarseth 2001). This threshold forinstability roughly agrees with the semi-major axis atwhich our binaries experience rapid eccentricity growthand eventual destruction. The slightly larger thresholdcompared to the theoretical expectation could be due tothe existence of disk, which leads to the onset of insta-bility more gently.3.5. BBH Accretion
Steady-state accretion
In steady-state the time and azimuthally averaged ac-cretion rate through the disk should be constant. Foran accreting BBH, the accretion rates in the disk inte-rior ( r < a com ) and exterior ( r > a com ) to the BBHCOM will be different but satisfy ˙ m in − ˙ m outer = ˙ m bin .Here a positive ˙ m bin denotes mass removal and a positive˙ m in , ˙ m outer refers to outwards mass transport. We findthat ˙ m bin depends on the adopted removal rate and thatif this rate is too fast the global flow of the disk can bereversed, i.e., ˙ m in can be outwards rather than inwards.We show this effect in Figure 5 by plotting the time-averaged ˙ m ( r ) = (cid:82) dφr Σ v r for two different viscositiesand domain sizes. The time-averaging is done over 100COM orbits at a time just before the BBH is released.For removal rates of τ − = Ω bh ,i and 5Ω bh ,i , ˙ m in >
0. In
Figure 5.
Disk accretion profiles for different removal ratesand disk viscosities. The blue line shows the azimuthal andtime averaged (between 400 −
500 orbits) disk accretion ratebefore the binary release for Model A in Table 1 as a functionof radial distance from the central SMBH (Note that thewave damping boundary layers have been cut off.). The blackline corresponds to Model Ar in Table 1. The jump at r = r matches exactly the accretion rate onto the BBH shown inFigure 7. The other two lines are with a higher disk viscosity( α = 1 . × − ) and two removal rates as indicated by thelegends. We can see that a large removal rate results in theinner disk accretion rate becomes unphysically positive. these situations, ˙ m bin > | ˙ m outer | , i.e the BBH wants toaccrete faster than the outer disk can feed it and so theinner disk must compensate. For a slower removal rateof 0 . bh ,i , the disk maintains inwards mass transporteverywhere and maintains ˙ m bin < | ˙ m outer | . In partic-ular, we find that just before release the BBH accretes ∼
50% of ˙ m outer for slow removal rates.3.5.2. Time evolution
While the global accretion rate in the disk sensitivelydepends on the removal rate, the time evolution of theBBH itself does not. This is shown in Figure 6 whichplots the evolution of the BBH for different accretionprescriptions. Despite the vastly different ˙ m ( r ) and˙ m bin , at early times the BBH expands and gains eccen-tricity at roughly the same rate – indicating that theselarge scale effects are unimportant near the BBH. Infact, we find that the flow and spatial torques on theCSD scale are roughly the same. If we were to takeinto account, the additional mass and angular momen-tum accreted by the binary, this result may change, butsimulations of isolated binaries find this component ofthe torque is typically small (Mu˜noz et al. 2019; Moodyet al. 2019).Figure 7 shows the time evolution of ˙ m bin and the ac-cretion rates of each individual BH for the slow removal Li et al.
Figure 6.
Same as Figure 2, but for different accretionmethods. Shown are Models A (black), Ar (magenta), andD (green). All three models show similar binary dynamics. case. At early times, the binary accretes symmetricallywith ˙ m bin ∼ [0 . − . | ˙ m outer | , with a periodic varia-tion of 2Ω bin , where Ω bin is the current binary orbitalfrequency.At later times when the binary has become eccentric,we find that one component of the binary preferentiallyaccretes more material than the other on a ∼ bin / t = 1250 COM orbits, thereis a noticeable slowdown in the BBH expansion rate for Figure 7.
BBH accretion rates for Model Ar. Upper panel:time-averaged BBH accretion rates. Red and blue lines cor-respond to accretion rates for two members of the binary,and the black line is their sum. The green line shows thetime evolution of the supplied accretion rate from the outerdisk (i.e., − ˙ m outer at r = 1 . r ). The lower panel showsa zoom-in of the accretion rates without averaging between1300 and 1301 orbits when the binary eccentricity is ≈ . our slow accretion model. We suspect that this couldbe related to the rapid increase of binary eccentricity,as simulations of isolated binaries have found that theexpansion rates for eccentric binaries are slower thantheir circular counterparts (Mu˜noz et al. 2019).3.6. Other Models: Influence of Binary InitialSeparation and Gap Depth
Initial separation — We have another two runs with theinitial separation of the binary of a = 0 . r (ModelsE, F). The binary evolution profiles are shown in Fig-ure 8. Similar to the Models A & B shown in Figure 2,two softenings with (cid:15) = 0 . a and (cid:15) = 0 . a have beentested. We obtain very similar results as the a = 0 . r cases: the binary is outspiraling when (cid:15) = 0 . a andit is inspiraling when (cid:15) = 0 . a , as long as the ini-tial separation is less than the critical value required totrigger the instability for the triple system discussed inSection 3.4. Gap depth and R H — We have tested several other caseswith a smaller gap depth controlled by K = (2 q ) /αh (e.g., two high α = 1 . × − models) and different bi-nary embeddedness parameters R H /H (e.g., R H /H (cid:46) R H /H > . − . a ), we always find that the binary expands BH in AGN Disks Figure 8.
Binary separation (shaded regions) and semi-major axis (solid lines) evolution for models with an initialseparation of a = 0 . r . Green ( (cid:15) = 0 . a , Model E)and blue ( (cid:15) = 0 . a , Model F) lines are for prograde orbitsand black ( (cid:15) = 0 . a , Model G) and red ( (cid:15) = 0 . a , ModelH) lines are retrograde orbits. Note that for the retrogradeorbits both small and large softening cases cause the binaryorbits contract with time. with time after the release, regardless of gap depth or R H /H . We have also tested the case of a similar small H as Tiede et al. (2020), which still shows an outspi-raling BBH, different from their inspiraling results. Athorough parameter study will be performed in the fu-ture. 3.7. Retrograde Binary
Up until now, we have only focused on prograde bina-ries. Here we further examine the impact of retrogradeorbits, i.e, when the BH-BH angular momentum is mu-tually inclined to that of the BBH-COM by 180 ◦ . Theseare Models G and H in Table 1.We show a close-up of the gas surface density andvelocity for Model G at 500 orbits in Figure 9. Theglobal gas surface density is quite similar to the progradecase shown in Figure 1. However, the gas distributionand flow pattern in the vicinity of the BBH are quitecomplex. There are clear CSDs that are aligned with theBBH’s angular momentum vector, but the spiral patternis much less pronounced than in the prograde case.We show the semi-major axis evolution for Models Gand H in Figure 8. Interestingly, the binary is inspi-raling even for a small softening. The inspiraling ratefor the retrograde orbit is larger than that of the largesoftening, prograde case (Baruteau et al. 2011). Fromthe separation evolution, the binary remains circular asit shrinks. Figure 9.
Same as Figure 1 but for Model G. The globalgas surface density is similar to that in Figure 1.
As we did for the prograde case, we show the cumu-lative torque and stacked, 2D torque density in Figures10 and 11. The surface density and torque density areaveraged over 400 snapshots between 500 − P com .Note that because the angular momentum vector of thebinary is flipped, we flip the sign of the torque so thatpositive values still correspond to increasing a bin .From Figure 10, we see that retrograde BBHs inspiralregardless of softening. In contrast to the prograde case,the CSD region actually makes the binary inspiral at afaster rate. By looking at Figure 11, we see that thisis because the CSDs lack any coherent spiral structurethat can provide an outspiraling torque. CONCLUSIONS AND DISCUSSIONWe have performed a series of global 2D isothermalhydrodynamical simulations of an equal mass-ratio BBHembedded in an accretion disk to study the evolution ofthe binary separation.We have shown that prograde binaries, which are mostoften discussed in the literature, expand rather than con-0
Li et al.
Figure 10.
Same as Figure 3 but for the retrograde orbitswith an initial separation of 0 . r for the binary. Dashedline indicates the CSD region with δr ≤ a . Note that − dT g ,<δr is shown here for the retrograde cases. tract, as is often assumed. This is because when theCSD region is adequately resolved – with both high res-olution and small softening – the positive CSD contri-bution to the total torque overwhelms the negative con-tribution from the CBD and inter-spiral regions. There-fore, a small softening with a high resolution to prop-erly resolve the CSD region is crucial, in agreement withrecent isolated binary simulations (Moody et al. 2019;Mu˜noz et al. 2019, 2020; Duffell et al. 2020; Tiede et al.2020). Such a conclusion does not sensitively dependon the accretion prescription for the binary we adopt.These results caution against the popularity of the AGNdisk channel for prograde BBH mergers as suggested inmany previous works (McKernan et al. 2012, 2014, 2018,2019; Bartos et al. 2017b; Stone et al. 2017; Leigh et al.2018; Secunda et al. 2019, 2020; Yang et al. 2019a,b;Gr¨obner et al. 2020; Ishibashi & Gr¨obner 2020; Abbottet al. 2020b; Tagawa et al. 2020a,b).However, when the softening scale is about half of theinitial binary separation (i.e., (cid:15) = 0 . a ), the binary or-bit contracts with time and stalls at the softening scale,consistent with Baruteau et al. (2011). The stalling ofthe orbital hardening at this large softening scale (cid:15) thuscannot bring the binary into a gravitational-wave (GW)-emitting dominated regime.In contrast to the prograde case, we find that retro-grade binaries do contract, regardless of softening. Thisis due primarily to the loss of the spiral structure inthe CSDs that provides the main source of expansion.This suggests that the gas-driven inspiral for retrograde binaries may produce a population of compact BBHsin the GW-emitting regime in AGN disks, which maycontribute a large fraction to the observed LIGO/Virgoevents. (as suggested by e.g., Secunda et al. 2020).For the outspiraling case, the binary’s eccentricitycan be significantly excited within a few hundred bi-nary COM orbits. Due to the eccentricity growth, theminimum separation between the binary could still besmaller than its initial value, both of which can enhancethe GW inspiraling rate significantly due to the strongdependence of the coalescence time scale on the binary’seccentricity (i.e., t gw ∝ (1 − e bin ) . ). We expect that thethe maximum binary’s eccentricity could be even higherif we decrease the initial semi-major axis of the binary.If this is the case, a highly eccentric BBH merger couldbe observed by future GW observations (Gayathri et al.2020; Samsing et al. 2020). It is also possible that thebinary will be destroyed due to the dynamical instabil-ity of the triple system, if they have not entered intothe GW-emitting regime yet after the rapid eccentricitygrowth stage.4.1. Limitations and Future Prospect
In this work we assume a small gravitational soft-ening by analogy with the isolated binary simulations.However, simulations of single intermediate mass ratiosecondaries use a much larger softening of (cid:15) (cid:39) . H .This number has been shown to produce approximatelyequal 2D and 3D gravitational torques (M¨uller et al.2012). However, for a 2D BBH simulation, this soften-ing is more appropriate for scales larger than a bin , i.e.,to obtain an equivalent 3D torque on the binary’s COMoutside of δr ∼ H , but certainly not for smaller scaleswhere the dynamics are very different from that of a sin-gle object. A study similar to M¨uller et al. (2012) whichcompares the 2D and 3D results for embedded BBHs isthus required to determine the “correct” softening scale.We suspect that the answer will actually be two soften-ing lengths, one for the larger scale that is related tothe COM potential, and one for the smaller scale thatis related to the higher moments of the BBH potential.Further work, however, is required.In addition to the question of softening, 3D simula-tions can also more accurately capture the vertical flowpattern around the CSD region and influence how thefully embedded binary evolves due to the 3D effect. Thedisk breaking due to the inclined binary with respect tothe embedded disk may also play some roles for the bi-nary dynamics.In all of our simulations, we have adopted a locallyisothermal equation of state, which could be unreal-istic, especially within the region of each BH’s Bondi BH in AGN Disks Figure 11.
Same as Figure 4 but for Models G & H which have retrograde orbits. The upper panels are for the softening of (cid:15) = 0 . a (Model G), the lower panels are for (cid:15) = 0 . a (Model H). Note that a minus sign is added to the torque density inthe middle and right panels for the retrograde cases. sphere. We intend to implement a more detailed treat-ment of the thermodynamics and/or more a realisticviscosity/temperature prescription around each BH totake account into the radiative cooling/heating of thegas. Souza Lima et al. (2017); del Valle & Volonteri(2018) highlighted that the AGN feedback can signif-icantly decrease the shrinkage rate of the SMBHB. Amore appropriate treatment of AGN radiative and me-chanical feedback physics (Yuan et al. 2018; Yoon et al.2018; Li et al. 2018) could better quantify their effect inthe binary evolution.Finally, there are various parameter spaces to be ex-plored in the future. This includes the binary mass ra-tio, the mass ratio to the central SMBH, the disk as-pect ratio, and disk mass. We should expect the binaryto inpiral for low mass ratio based on the consensus ofinward migration for a giant planet in protoplanetarydisks. The transition from outspiraling to inspiralingfor a decreasing mass ratio in isolated binary simulationhas been reported recently (Duffell et al. 2020, Dempseyet al. 2020, in prep.). This will lead to an extreme massratio inspiral event from the BBH merger. We have explored a regime with the disk mass being less thanthe mass of BBH, but the systems with more massivedisk could behave very differently (Cuadra et al. 2009;Roedig et al. 2012), which may also drive the binary tocoalescence. ACKNOWLEDGMENTSWe would like to thank Douglas Lin, Cl´ementBaruteau, Bin Liu for beneficial discussions. A.M.Dthanks Diego Mu˜noz for many helpful discussions onbinary-disk interaction. We gratefully acknowledge thesupport by LANL/LDRD and NASA/ATP. This re-search used resources provided by the Los Alamos Na-tional Laboratory Institutional Computing Program,which is supported by the U.S. Department of EnergyNational Nuclear Security Administration under Con-tract No. 89233218CNA000001. The LANL Reportnumber is LA-UR-20-28080. Softwares: LA-COMPASS (Liet al. 2005, 2009),
Numpy (van der Walt et al. 2011),
Matplotlib (Hunter 2007)2
Li et al.
REFERENCES
Abbott, B. P., Abbott, R., Abbott, T. D., LIGO ScientificCollaboration, & Virgo Collaboration. 2016, PhRvL, 116,061102Abbott, B. P., Abbott, R., Abbott, T. D. e. a., LIGOScientific Collaboration, & Virgo Collaboration. 2019,Physical Review X, 9, 031040Abbott, R., Abbott, T. D., Abraham, S., & Acernese, F.e. a. 2020a, arXiv e-prints, arXiv:2010.14527Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b,ApJL, 900, L13Antonini, F., Toonen, S., & Hamers, A. S. 2017, ApJ, 841,77Banerjee, S. 2017, MNRAS, 467, 524Bartos, I., Haiman, Z., Marka, Z., et al. 2017a, NatureCommunications, 8, 831Bartos, I., Kocsis, B., Haiman, Z., & M´arka, S. 2017b, ApJ,835, 165Baruteau, C., Cuadra, J., & Lin, D. N. C. 2011, ApJ, 726,28Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ,714, 1217Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R.2016, Nature, 534, 512Bellovary, J. M., Mac Low, M.-M., McKernan, B., & Ford,K. E. S. 2016, ApJL, 819, L17Chen, X., Lin, D. N. C., & Zhang, X. 2020a, ApJL, 893, L15Chen, Y.-X., Zhang, X., Li, Y.-P., Li, H., & Lin, D. N. C.2020b, ApJ, 900, 44Cresswell, P., & Nelson, R. P. 2006, A&A, 450, 833Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181,587Cuadra, J., Armitage, P. J., Alexander, R. D., & Begelman,M. C. 2009, MNRAS, 393, 1423de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al.2006, MNRAS, 370, 529del Valle, L., & Volonteri, M. 2018, MNRAS, 480, 439Dempsey, A. M., Lee, W.-K., & Lithwick, Y. 2020, ApJ,891, 108Derdzinski, A., D’Orazio, D., Duffell, P., Haiman, Z., &MacFadyen, A. 2021, MNRAS, 501, 3540D’Orazio, D. J., Haiman, Z., & MacFadyen, A. 2013,MNRAS, 436, 2997Duffell, P. C., D’Orazio, D., Derdzinski, A., et al. 2020,ApJ, 901, 25Dunhill, A. C., Cuadra, J., & Dougados, C. 2015, MNRAS,448, 3545Eggleton, P., & Kiseleva, L. 1995, ApJ, 455, 640Farris, B. D., Duffell, P., MacFadyen, A. I., & Haiman, Z.2014, ApJ, 783, 134 Fragione, G., & Kocsis, B. 2019, MNRAS, 486, 4781Gayathri, V., Healy, J., Lange, J., et al. 2020, arXive-prints, arXiv:2009.05461Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425Graham, M. J., Ford, K. E. S., McKernan, B., et al. 2020,PhRvL, 124, 251102Gr¨obner, M., Ishibashi, W., Tiwari, S., Haney, M., &Jetzer, P. 2020, A&A, 638, A119Heath, R. M., & Nixon, C. J. 2020, A&A, 641, A64Hunter, J. D. 2007, Computing in Science and Engineering,9, 90Ishibashi, W., & Gr¨obner, M. 2020, A&A, 639, A108Kim, H., & Kim, W.-T. 2007, ApJ, 665, 432Kim, H., Kim, W.-T., & S´anchez-Salcedo, F. J. 2008,ApJL, 679, L33Leigh, N. W. C., Geller, A. M., McKernan, B., et al. 2018,MNRAS, 474, 5672Li, H., Li, S., Koller, J., et al. 2005, ApJ, 624, 1003Li, H., Lubow, S. H., Li, S., & Lin, D. N. C. 2009, ApJL,690, L52Li, Y.-P., Chen, Y.-X., Lin, D. N. C., & Zhang, X. 2021,ApJ, 906, 52Li, Y.-P., Yuan, F., Mo, H., et al. 2018, ApJ, 866, 70Lin, D. N. C., & Papaloizou, J. C. B. 1993, in Protostarsand Planets III, ed. E. H. Levy & J. I. Lunine, 749Liu, B., & Lai, D. 2017, ApJL, 846, L11—. 2019, MNRAS, 483, 4060MacFadyen, A. I., & Milosavljevi´c, M. 2008, ApJ, 672, 83Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634Mardling, R. A., & Aarseth, S. J. 2001, MNRAS, 321, 398McKernan, B., Ford, K. E. S., Kocsis, B., Lyra, W., &Winter, L. M. 2014, MNRAS, 441, 900McKernan, B., Ford, K. E. S., Lyra, W., & Perets, H. B.2012, MNRAS, 425, 460McKernan, B., Ford, K. E. S., Bellovary, J., et al. 2018,ApJ, 866, 66McKernan, B., Ford, K. E. S., Bartos, I., et al. 2019, ApJL,884, L50Michaely, E., & Perets, H. B. 2019, ApJL, 887, L36Miranda, R., Mu˜noz, D. J., & Lai, D. 2017, MNRAS, 466,1170Moody, M. S. L., Shi, J.-M., & Stone, J. M. 2019, ApJ,875, 66Mu˜noz, D. J., & Lai, D. 2016, ApJ, 827, 43Mu˜noz, D. J., Lai, D., Kratter, K., & Mirand a, R. 2020,ApJ, 889, 114Mu˜noz, D. J., Miranda, R., & Lai, D. 2019, ApJ, 871, 84M¨uller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541,A123
BH in AGN Disks13