Eccentric Black Hole Mergers in Active Galactic Nuclei

The astrophysical origin of gravitational wave transients is a timely open question in the wake of discoveries by the Laser Interferometer Gravitational-Wave Observatory (LIGO)/Virgo. In active galactic nuclei (AGNs), binaries form and evolve efficiently by interaction with a dense population of stars and the gaseous AGN disk. Previous studies have shown that stellar-mass black hole (BH) mergers in such environments can explain the merger rate and the number of suspected hierarchical mergers observed by LIGO/Virgo. The binary eccentricity distribution can provide further information to distinguish between astrophysical models. Here we derive the eccentricity distribution of BH mergers in AGN disks. We find that eccentricity is mainly due to binary–single (BS) interactions, which lead to most BH mergers in AGN disks having a significant eccentricity at 0.01 Hz, detectable by the Laser Interferometer Space Antenna. If BS interactions occur in isotropic-3D directions, then 8%–30% of the mergers in AGN disks will have eccentricities at 10 Hz above e10 Hz ≳ 0.03, detectable by LIGO/Virgo/Kamioka Gravitational Wave Detector, while 5%–17% of mergers have e10 Hz ≥ 0.3. On the other hand, if BS interactions are confined to the AGN–disk plane due to torques from the disk, with 1–20 intermediate binary states during each interaction, or if BHs can migrate to ≲ 10−3 pc from the central supermassive BH, then 10%–70% of the mergers will be highly eccentric (e10 Hz ≥ 0.3), consistent with the possible high eccentricity in GW190521.


INTRODUCTION
Recent detections of gravitational waves (GWs) have shown evidence for a high black hole (BH) merger rate in the Universe (Abbott et al. 2019a).However, the proposed astrophysical pathways to merger remain debated.
Recently, the merger of two unusually heavy BHs, GW190521, was reported (Abbott et al. 2020a,b).Their high masses (85 +21 −21 M and 66 +17 −18 M ) and high spins (0.69 +0.27  −0.62 and 0.73 +0.24  −0.64 ) indicate that the merging BHs themselves could have been the remnants of earlier BH E-mail: htagawa@astr.tohoku.ac.jp mergers, which increased their masses beyond the ∼ 56 M limit due to pulsational pair-instability (Farmer et al. 2019) and their spins beyond the small natal values expected in stellar evolutionary models of angular momentum transfer (Fuller & Ma 2019).Thus, detection of GW190521 (and also other events: GW170729, GW170817A, GW190412 and GW190814; Zackay et al. 2019;Abbott et al. 2020c,d) suggests that hierarchical mergers could be frequent among compact objects.This is consistent with the scenario of mergers in AGN disks (e.g.Yang et al. 2019Yang et al. , 2020;;Tagawa et al. 2020a, hereafter Paper II), in which binaries are efficiently hardened by interaction with the surrounding gas (e.g.Bartos et al. 2017;Stone et al. 2017;McKernan et al. 2018) and with the dense populations of stars and compact objects (Tagawa et al. 2020b, hereafter Paper I).A possible electromagnetic counterpart would further support this scenario (McKernan et al. 2019;Graham et al. 2020).
Interestingly, Gayathri et al. (2020) found that GW190521 prefers a high eccentricity of e 10Hz ∼ 0.7, along with high spin-precession (see also Romero-Shaw et al. 2020).A high eccentricity places additional constraints on possible astrophysical models of this source (e.g.Zevin et al. 2019;Rodriguez et al. 2018;Rasskazov & Kocsis 2019).Samsing et al. (2020b) recently showed that highlyeccentric mergers are common in AGN disks by assuming that binary-single (BS) interactions of BHs are confined to a plane.That study assumed constant values for the frequency of BS interactions, initial separation, and relative velocity of a third body based on Paper I.
In this Letter, we investigate the distribution of the binary eccentricity for mergers in AGN disks by performing one dimensional (1D) N -body simulations combined with semi-analytical prescriptions, which enable us to follow binaries considering such effects as eccentricity evolution due to BS interactions, type I/II torques exerted by circumbinary disks, and GW radiation ( § 2.3).We also augment the model used in Papers I/II to include GW capture in single-single (SS) and BS encounters.
2. METHOD Our model is based on that in Paper I and partially in Paper II as specified below.The new ingredients are described in some detail in § 2.2.1, § 2.3.2 and § 2.4.

Components
We suppose that there is a supermassive BH (SMBH) at the center of a galaxy and it has a gaseous accretion disk (hereafter "AGN disk") and a spherically distributed stellar cluster (hereafter "spherical component").We follow the evolution of the BH system consisting of a cluster of BHs around the AGN disk (called the flattened component due to its spatial distribution, which is assumed to be caused by vector resonant relaxation; Szolgyen & Kocsis 2018), and stars and BHs captured inside the AGN disk due to the gaseous torque (called the disk stellar/BH components).See Fig. 1 in Paper I for an illustration of these components.

Binary formation and disruption
Some BHs in the flattened component are in binaries from the outset (called pre-existing binaries).In the AGN disk, binaries form due to gas dynamical friction during two-body encounters (dubbed the gas-capture binary formation).Binaries also form due to dynamical interactions during three-body encounters (dynamical binary formation).At their formation, a thermal distribution f (e) = 2e is assumed, while e = 0 was assumed in Papers I/II.In addition to those mechanisms already included in Paper I, we here consider the binaries formed by the GWC mechanism (e.g.O' Leary et al. 2009;Samsing et al. 2020a), which are relevant for highly eccentric mergers, as described below in § 2.2.1.Binaries are disrupted by soft-BS interactions, in which the binaries become softer (e.g.Heggie 1975).

Binary formation by the GW capture
We treat binary formation by GWC due to SS encounters (SS-GWC) and BS interactions (BS-GWC) separately as described below.
In an SS encounter event, a binary can form if the two bodies approach close enough that the energy radiated by GWs (∆E GW ) exceeds their kinetic energy, E SS ≈ 1 2 µv 2 12 , where µ is the reduced mass, and v 12 the relative velocity.∆E GW is approximately given as (e.g.O' Leary et al. 2009), where G is the gravitational constant, c the light speed, m tot the total mass of the two bodies, η the symmetric mass ratio, and the pericenter distance where b is the impact parameter of the encounter.By equating ∆E GW in equation ( 1) and E SS , the maximum pericenter distance for the GWC is Note that the gas has negligible impact on the dynamics during the GWC owing to very short timescale for the capture (∼ hr, O' Leary et al. 2009).Here we define the maximum impact parameter b max at which r p = r p,max .We assume that the cross section of SS-GWC σ SS is approximately given as σ SS = r SS z SS , where r SS = min[b max , r Hill ], z SS = min[r SS , h c ], r Hill = r(m BH /3M SMBH ) 1/3 is the Hill radius with respect to the SMBH, m BH is the mass of the BH, and h c is the average orbital height of background objects at the radial location of the BH (see Paper I).
For each object, the timescale to undergo the SS-GWC is t GSS = 1/(n c σ SS v rel ), where n c is the number density of background objects, and v rel the typical relative velocity between the BH and background objects (Paper I).Here we substitute v rel = v 12 in Eqs. ( 1)-(2).We set the probability of binary formation by this mechanism during the timestep ∆t to P GSS = ∆t/t GSS .When the SS-GWC binaries form, we choose b so that b 2 is uniformly distributed between 0 and b 2 max .The semi-major axis and e of the newly formed binary are s = Gm 2 tot η/(2|∆E GW | − 2E SS ), and e = 1 − r p /s, respectively.We only consider binary formation through SS-GWC between disk BHs, as the probability of interactions is much lower for spherical component objects due to their large v rel (P SS ∝ v −11/7 rel for b max < h c , r Hill ), and due to the large physical sizes of stars.For binaries formed in this process, we assume that the orbital angular momentum directions are aligned (or half antialigned) with the angular momentum direction of the AGN disk.
GW-capture binaries can also form during BS interactions (e.g.Samsing et al. 2014) if one body is captured due to strong GW emission.Before the interaction, the velocities among BHs in the flattened component are likely distributed almost randomly in 3D space, while for interactions among the disk BHs the three-body velocities may be constrained in the AGN-disk plane due to the disk gas drag.We here study the 3D and 2D hard-BS interactions separately.In the case of the 3D hard-BS interactions, the probability of GW-capture binary formation is roughly given by (Samsing 2018), where N int is the average number of quasi-stable states for temporary binary BHs formed during a BS interaction, where N int ∼ 20 for equal-mass cases (Samsing et al. 2014;Samsing et al. 2020b).The GWC probability is calculated as  is the probability that GWC occurs in one interaction, e cr = 1 − r p,max /s is the critical eccentricity above which the GWC binaries form due to GW emission, and f (e) is either f 3D = 2e (e.g.Heggie 1975) or Valtonen & Karttunen 2006) for isotropic-3D and 2D interactions, respectively.
When the BS-GW-capture binaries form, assuming the semi-major axis s does not change from that before the BS interaction, we determine r p = s(1 − e) by drawing e from f (e).Capture occurs only if r p < r p,max (Eq.3).We assume that v 12 = Gm tot /s, i.e. the intrinsic orbital velocity of the binary.For simplicity, we also assume that the BS-GWC results in binary formation among the most and second-most massive objects.Although a BS-GW-capture binary may be bound to the third object if the GW kick is small (Samsing & Ilan 2019), we ignore the third body after the binary formation for simplicity.

Gaseous processes
The velocity of BHs relative to the local motion of the AGN disk decreases due to the accretion torque and gas dynamical friction.The semi-major axis (s) evolves due to gas dynamical friction by the AGN disk and type I/II migration torques exerted by a compact circumbinary disk that forms within the Hill sphere of the binary.The radial distances of BHs from the SMBH (r) are also allowed to evolve due to type I/II torques from the AGN disk.Gas accretion affects BH masses, BH spins, and the orbital angular momentum directions of binaries (Paper II).When gaps form around BHs, we assume that migration torques, gas dynamical friction, and gas accretion are weakened.Objects usually accumulate in gap forming regions (Paper I), which act like migration traps.

Dynamical processes
We also account for dynamical interactions with single stars/BHs and BH binaries in a way similar to the Monte-Carlo method (Paper I).The binaries' semi-major axes, velocities, and orbital angular momentum directions evolve due to BS interactions, and the velocities of all BHs additionally evolve due to scattering.When the interactions are confined in 2D, kicks are also assigned within the AGN-disk plane.
While binary component exchanges were neglected during BS interactions in Papers I/II, here we assume that the components are always replaced by the most massive pair during hard-BS interactions, in which the binary becomes harder.During an exchange, the binary's binding energy is assumed to be unchanged (Sigurdsson & Phinney 1993).See the Appendix for the impact of including these exchange interactions.In binary-binary interactions, we assume that the binary composed of the two most massive objects experiences two BS interactions with the two less massive objects, which become unbound after the interaction.

Merger prescription
Once binaries become sufficiently tight, their separation shrinks promptly by GW emission.We assume that the BHs merge when their pericenter distance becomes smaller than the innermost stable circular orbit, and then assign a kick velocity and mass loss due to the GW emission and prescribe BH spin evolution at the merger as in Paper II.

Binary eccentricity evolution
After formation, the binary eccentricity e changes due to GW emission, the torque by the circumbinary disk, and BS interactions.Below we describe our prescriptions for calculating those processes, which are newly incorporated in our model.
For the GW emission, we follow Peters (1964) to track the change of the eccentricity.
For the torque by the circumbinary disk, we use the fitting function of Eq. ( 2) of Zrake et al. ( 2020) in e ≤ 0.8 and assume de/dlogm bin = −2.3 in e > 0.8.Their results using 2D hydrodynamical simulations suggest disks drive binaries to e ≈ 0.45 (see also Muñoz et al. 2019).
During hard-BS interactions, we draw e randomly from f (e) (see below Eq. 6).In the fiducial model, we assume that e follows f 3D (e) after all hard-BS interactions, but see different choices in the Appendix.

Models examined
The fiducial model parameters (M1) are listed in Table 1, which are the same as those in Papers I/II.We also study the additional models M2-M4, which differ in the treatment of the BS interactions.While 3D BS interactions are assumed in model M1, motions of the three bodies may be confined to the AGN-disk plane before the interaction.Even in the latter case, the motion may be isotropized in 3D during the interaction by small perturbations, such as inhomogeneities/warps of the AGN disk and presence of other objects.To mimic such 2→3D interactions, we examine another model (M2) in which the capture probability during a BS interaction is given by Eq. ( 5) assuming that randomization to 3D occurs after N int states, while the e distribution after the interaction obeys f 3D (e) and the binary angular momentum direction is randomized unless captured.The reason to adopt the model with 2→3D BS interactions instead of 2D BS interactions is stated in Appendix.To assess the importance of BS interactions, we omit it altogether in model M3.Similar to model M3, but more realistic is model M4, in which BS interactions are made inefficient by reducing the initial BH number (N BH,ini = 6×10 3 ) and enhancing the initial BH velocity dispersion (β v = 1).

RESULTS
In Fig. 1, we show the eccentricity distribution of mergers and the contribution of individual binary formation channels to it.Fig. 2 shows the expected distribution of the mass-weighted sum of the individual spin components perpendicular to the orbital plane χ eff , the precession parameter χ p (upper panels), the primary mass m 1 , and e (lower panels) of the observed mergers except for model M3.For M3, the χ eff and χ p distributions are not shown since χ p is always 0 due to the assumption that the orbital planes of SS-GWC binaries are aligned with the AGN disk plane ( § 2.2.1).The merger rate is weighted by the detection volume in the same manner as in Papers I/II and the eccentricity is calculated at the GW frequency of 10 Hz as in Eq. (36) of Wen (2003).
In the fiducial model M1, where the BS interaction is assumed to be 3D, ∼ 14% of mergers have a significantly nonzero eccentricity (e 10Hz 0.03) measurable by the current ground-based detectors LIGO/Virgo/KAGRA (e.g.Romero-Shaw et al. 2019, but see also Huerta & Brown 2013), due to GWC binary formation (Fig. 1 and Table 2) and frequent BS interactions.While most AGN merger events take place among binaries of non-GWC origin, which have lower e (Fig. 1 a), LISA observations will be able to measure e as low as e 10 −3 − 10 −2 at ∼ 0.01 Hz (e.g.Nishizawa et al. 2016), corresponding to e 10Hz 10 −7 if binaries evolve only by GW emission.e 10Hz are typically distributed in higher values (∼ 10 −3 -10 −2 ) than those of mergers in globular clusters (∼ 10 −7 -10 −4 ; Rodriguez et al. 2018) and isolated binaries (∼ 10 −7 -10 −6 ; Breivik et al. 2016), while the similar values to those in triple system (∼ 10 −3 -10 −2 ; Antonini et al. 2017).Due to randomization of the orbital angular momentum direction by frequent BS interactions before mergers, χ eff , and χ p tend to be distributed at low and high values (upper panel, Fig. 2 a), respectively.The physical properties, χ eff , χ p , and m 1 of GW190521 are generally consistent with the predictions by model M1, while the possible high e implied by Gayathri et al. (2020) is in tension with this model (lower panel, Fig. 2 a).
In model M2 with 2→3D BS interactions, the fraction of the BS-GWC binary mergers among the detectable events f GWBS is elevated to 0.84 (from 0.10 of model M1; see Table 2), and they have very high eccentricity e 0.1 (Fig. 1 b).In this way, BS interactions significantly increase the rate of highly eccentric mergers (see also Samsing et al. 2020b).The mass of GW190521, however, appears to be somewhat too high compared with the typical values in model M2 (lower panel b of Fig. 2), but the model may still be consistent with the event as the mass distribution for high-e events has been poorly constrained.
In model M3, where the BS interaction is omitted, the (SS-)GWC binary mergers occupy a high fraction of f GWSS = 0.54 favoring mergers between high-generation BHs (upper panel c of Fig. 2).High-generation BHs are those that have already experienced mergers in the past, and they tend to have high eccentricities around 0.1 (Fig. 1 c).Here, binary formation occurs mostly by the (SS-)GWC because BHs rapidly migrate without delay by kicks at the BS interactions to the inner AGN disk of r 10 −3 pc (e.g.Fig. 1 of Paper II), where non-GW-capture channels by the gas-capture and dynamical binary formation are inefficient due to high shear velocity ( tens km/s).
In general, the frequency of the SS-GWC f GWSS increases as the initial BH number decreases (see Appendix).For example, in model M4, where N ini = 6000 and β v = 1, f GWSS becomes as high as 0.38.The prediction of positive χ p , high e, and high-BH mass m 1 in this model agrees with the observed properties of GW190521 (Fig. 2d).Although the number of mergers is smaller by Table 2 The assumptions and results in different models.The input columns show the model number and the differences with respect to the fiducial model."BS in 2→3D" represents the model in which BS interactions occur in 2→3D space ( § 2.5), "No BS" represents the model in which BS interactions do not operate.The output columns list the detection fraction of mergers among non GWC binaries (f noGWC ), GWC binaries formed in SS encounters (f GWSS ), and those in BS interactions (f GWBS ) with respect to all mergers, the detection fractions of mergers whose e 10 Hz are in the range 0-0.03, 0.03-0.3,0.3-0.9, and 0.9-1, respectively, and the number of mergers (Nmer).The cumulative (upper) and differential (lower) detection rate distributions of the eccentricity at 10Hz e 10Hz for models M1 (a), M2 (b) in which BS interactions occur in 2→3D space ( § 2.5), M3 (c) in which BS interaction does not operate, and M4 (d) in which the initial BH number is low (N ini = 6000) and the initial velocity dispersion of BHs is high (βv = 1).In addition to the sum of all binary mergers (black), separate contributions from binaries due to the non-GWC (green), SS-GWC (orange), and BS-GWC (red) are shown.
a factor of ∼ 10 compared to that in the fiducial model, it is still consistent with the rate of GW190521-like events (see § 5.5 of Paper I).
As seen above, the AGN origin can naturally reproduce the occurrence of GW190521-like mergers.In several models, the eccentricity distribution allows us to distinguish specific types of the GWC events (orange and red lines in Fig. 1) due to differences in the initial relative velocity involved (Eq. 3) and such high eccentricities are measurable by LIGO/Virgo/KAGRA.The consistency of other observables of GW190521 will be examined in a follow-up paper.
Our calculations with different sets of model parameters of the AGN disk, BH population, etc. show that the predicted eccentricity distribution is robust, i.e. it is not very sensitive to these assumptions (see Appendix) other than those related to BS interactions as presented above.

CONCLUSIONS
In this Letter we have investigated the eccentricity distribution of merging black holes (BHs) in active galactic nuclei (AGN) accretion disks.Our conclusions are summarized as follows: 1. Due to frequent binary-single interactions in gapforming regions of the AGN disk, eccentricities at 10 Hz are above e 10Hz ≥ 10 −4 for all stellar-mass BH mergers in AGNs in our models.LIGO together with LISA will be able to constrain the corresponding astrophysical models at high and low frequencies, respectively.

If binary-single interactions occur in 2D or 2→3D
directions due to torques from the AGN disk, the mergers among gravitational wave capture binaries formed during BS interactions become very frequent (∼ 15-90%).Also, if BHs can migrate to r 10 −3 pc due to low BH density in the AGN disk, those formed in single-single encounters are significantly enhanced (∼ 40-70%).In these cases, ∼ 10-70% of detectable mergers have e 10Hz ≥ 0.3, which may explain the possible extreme eccentricity implied for GW190521 by Gayathri et al. (2020).
We have neglected the evolution of e by soft-BS interactions, weak scattering, and gas dynamical friction, which presumably have minor effects on e as their contributions on the later evolution of binaries are negligibly small.However, hierarchical triples, neglected in this study, may significantly affect the e distribution.We also have neglected other processes including the formation and evolution of compact objects (and their binaries) other than BHs, and the possible presence of massive perturbers (Deme et al. 2020), which may affect various properties of mergers.Also, BS interactions for nearly 2D cases need to be elucidated using N -body simulations.These effects will be investigated in future work.
which BS interactions occur isotropic in 3D Number of temporary binary BHs formed during a BS interaction N int = 20 Initial BH spin magnitude |a| = 0 Angular momentum directions of circum-BH disks ĴCBHD = ĴAGN for single BHs, ĴCBHD = Ĵbin for BHs in binaries Ratio of viscosity responsible for warp propagation over that for transferring angular momentum ν 2 /ν 1 = 10 Alignment efficiency of the binary orbital angular momentum due to gas capture frot = 1 (Eq.14 in Paper II) Mass of the central SMBH M SMBH = 4 × 10 6 M Gas accretion rate at the outer radius of the simulation (5 pc) Ṁout = 0.1 ṀEdd with η = 0.1 Fraction of pre-existing binaries fpre = 0.15 Power-law exponent for the initial density profile for BHs γρ = 0 Initial velocity anisotropy parameter such that βvv kep (r) is the BH velocity dispersion βv = 0.2 Efficiency of angular momentum transport in the α-disk α SS = 0.1 Stellar mass within 3 pc M star,3pc = 10 7 M Stellar initial mass function slope δ IMF = 2.35 Angular momentum transfer parameter in the outer star forming regions m AM = 0.15 (Eq.C8 in Thompson et al. 2005) Accretion rate in Eddington units onto stellar-mass BHs with a radiative efficiency η = 0.1 Γ Edd,cir = 1 Numerical time-step parameter ηt = 0.1 Number of radial cells storing physical quantities N cell = 120 Maximum and minimum r for the initial BH distribution r in,BH = 10 −4 pc, r out,BH = 3 pc Initial number of BHs within 3 pc N BH,ini = 2 × 10 4 Figure1.The cumulative (upper) and differential (lower) detection rate distributions of the eccentricity at 10Hz e 10Hz for models M1 (a), M2 (b) in which BS interactions occur in 2→3D space ( § 2.5), M3 (c) in which BS interaction does not operate, and M4 (d) in which the initial BH number is low (N ini = 6000) and the initial velocity dispersion of BHs is high (βv = 1).In addition to the sum of all binary mergers (black), separate contributions from binaries due to the non-GWC (green), SS-GWC (orange), and BS-GWC (red) are shown.

Figure 2 .
Figure2.Detection rate distributions of several parameters for the same models as in Fig.1.Upper and lower panels show distributions of χ eff -χp and m 1 -log 10 e 10Hz , respectively, while in model M3, m 1 -log 10 e 10Hz are shown for all mergers (lower) and mergers among n-generation (n ≥ 2) BHs (upper).The error bars correspond to the 90 percentile credible intervals for GW190521(Abbott et al. 2020a), while we here roughly assume e 10Hz = 0.7 +0.1 −0.2 referring toGayathri et al. (2020).

Figure 3 .
Figure3.The cumulative (upper) and differential (lower) detection rate distributions for the eccentricity at 10 HZ among all merging binaries for models M1, M2, M5-M9, M12-M16, and M18-M22.Different colors and rows present the distributions for different variations of the models, as labeled on the upper legends and listed in Table2.

Table 1
Fiducial values of our model parameters.

Table 3
Same with Table3, but results for all models."BS in 2D" represents the model in which BS interactions occur in 2D."No BS" represents the model in which BS interactions do not operate."No gas torque on ecc", "No gas hardening", and "No gas migration" respectively represent the models in which the gaseous torques are neglected with respect to the evolution of e, the binary semi-major axis, and the radial position of the binary center of mass within the AGN disk."3 × m 1g " represents the model in which initial BH masses are multiplied by 3. "No exchange in BS" represents the model in which the binary components are not exchanged during BS interactions.