Spin Evolution of Stellar-mass Black Hole Binaries in Active Galactic Nuclei

The astrophysical origin of gravitational wave (GW) events is one of the most timely problems in the wake of the LIGO/Virgo discoveries. In active galactic nuclei (AGN), binaries form and evolve efficiently by dynamical interactions and gaseous dissipation. Previous studies have suggested that binary black hole (BBH) mergers in AGN disks can contribute significantly to BBH mergers observed by GW interferometers. Here we examine the distribution of the effective spin parameter $\chi_\mathrm{eff}$ of this GW source population. We extend our semi-analytical model of binary formation and evolution in AGN disks by following the evolution of the binary orbital angular momenta and black hole (BH) spins. BH spins change due to gas accretion and BH mergers, while the binary orbital angular momenta evolve due to gas accretion and binary-single interactions. We find that the distribution of $\chi_\mathrm{eff}$ predicted by our AGN model is similar to the distribution observed during LIGO/Virgo O1 and O2. On the other hand, if radial migration of BHs is inefficient, $\chi_\mathrm{eff}$ is skewed toward higher values compared with the observed distribution, because of the paucity of scattering events that would randomize spin directions relative to the orbital plane. We suggest that high binary masses and the positive correlation between binary mass and the standard deviation of $\chi_\mathrm{eff}$ for chirp masses up to $\approx 20\,\mathrm{M}_\odot$, can be possible signatures for mergers originating in AGN disks. Finally, hierarchical mergers in AGN disks naturally produce properties of the recent GW event GW190412, including a low mass ratio, a high primary BH spin, and a significant spin component in the orbital plane.

Galactic nuclei are the densest environments of stars E-mail: htagawa@caesar.elte.hu and compact objects in the Universe (see Neumayer et al. 2020, for a recent review). In an active galactic nucleus (AGN), a high-density gas disk forms within 0.1-10 pc (Burtscher et al. 2013) around a central supermassive BH (SMBH). Several authors have recently pointed out that these environments are conducive to forming compact-object binaries. This "AGN disk channel" has received increasing attention in the wake of the LIGO/Virgo discoveries, as a possible explanation for some of the LIGO/Virgo events. In particular, McKernan et al. (2012McKernan et al. ( , 2014 predicted the formation of intermediate-mass BHs in AGN disks due to collisions of compact objects. Bartos et al. (2017b) have proposed a pathway for binary BH (BBH) mergers in AGN disks in which binaries are captured by an accretion disk within ∼ 0.01 pc from the SMBH due to linear momentum exchange during disk-crossing, and after that, binaries are hardened by gas dynamical friction by an AGN disk and type I/II torque by circumbinary disks. Stone et al. (2017) have proposed another pathway, in which in-situ formed binaries at ∼pc scale evolve via effects of binarysingle interactions with a disk stellar component and type I/II torque from circumbinary disks. Bellovary et al. (2016) suggested that BHs accumulate and merge with each other in migration traps at 20 − 300 Schwarzschild radii where the sign of the torque from the AGN disk changes. Secunda et al. (2018), Yang et al. (2019a,b) and Gayathri et al. (2019) investigated the properties of mergers in migration traps. Tagawa et al. (2019) investigated how binaries form and merge in AGN disks by performing self-consistent one-dimensional N -body sim-ulations combined with semi-analytical prescriptions of the relevant processes. They found that binaries form efficiently in the inner regions ( pc; but outside the migration traps) of AGN disks, due to the dissipation of relative velocities of unbound pairs of BHs via gas drag ("gas-capture" binaries). Motivated by the above, in the present paper, we investigate whether the AGN disk channel can be distinguished from other formation pathways. Previous work proposed distinguishing features based on spatial associations with bright AGN (Bartos et al. 2017a;Corley et al. 2019), large chirp masses (Tagawa et al. 2019), redshift evolution (Yang et al. 2020), acceleration of the binary's center of mass (Meiron et al. 2017;Inayoshi et al. 2017b;Wong et al. 2019), or gravitational lensing (Kocsis 2013;D'Orazio & Loeb 2019). One feature that has not yet been studied in detail in this channel is the expected distribution of BH spins. The effective spin parameter χ eff (which is the sum of the projection of binary spins onto the binary orbital angular momentum) has been shown to provide useful information to constrain other compact-object merger pathways Vitale et al. 2017;Talbot & Thrane 2017).
The χ eff distribution inferred from the observed GW events prefers low values (The LIGO Scientific Collaboration et al. 2018), which suggests low natal BH spins or random directions between the binary orbital angular momentum and the BH spins (Farr et al. 2017). On the other hand, several events are reported to have high or low χ eff values (Zackay et al. 2019b,a; The LIGO Scientific Collaboration & the Virgo Collaboration 2020). Safarzadeh et al. (2020) suggested that there is a negative and positive correlation between mass and the mean and the dispersion of χ eff , respectively. For the evolution of isolated binaries, the low observed χ eff values may be reproduced if the angular momentum transport within the stars is highly efficient (Qin et al. 2018;Bavera et al. 2019). In globular clusters, the orbital angular momentum directions of binaries are randomized by binary-single interactions, which predicts a χ eff distribution symmetric around χ eff = 0 and favoring low values (Rodriguez et al. 2018;Arca-Sedda et al. 2018). For triple systems, the Kozai mechanism can cause BH spin misalignment (Liu & Lai 2017Liu et al. 2019;Fragione & Kocsis 2019).
For mergers in AGN disks, Yang et al. (2019b), McKernan et al. (2019), and Gayathri et al. (2019) estimated the χ eff distribution for mergers in migration traps, and McKernan et al. (2019) estimated the spin distribution for mergers among binaries formed during close encounters in AGN disks, by assuming that binaries are always aligned or anti-aligned with the AGN disk, and BH mergers are much faster than the growth of BH spins by gas accretion. However, Tagawa et al. (2019, hereafter Paper I) have suggested that the orbital angular momenta of the binaries at mergers are often misaligned with the AGN disk due to frequent hard binary-single interactions. Since the merger timescale of BBHs is comparable to the timescale of gas accretion onto BHs, the evolution of the BH spins needs to be explicitly followed, accounting for both effects. Furthermore, to determine the χ eff distribution for mergers in AGN disks, it is necessary to also follow the orbital angular momenta of the binaries, again taking into account both binary-single interactions and gas accretion.
In this paper, we determine the distribution of χ eff for binary mergers in AGN disks, by incorporating the evolution of BH spins and the binary orbital angular momenta into the semi-analytical prescriptions and onedimensional N-body simulations used in Paper I. We find that the frequency of binary-single interactions, the angular momentum of the captured gas, and the initial BH spin directions strongly influence the χ eff distribution. The rest of this paper is organized as follows. In § 2, we describe the numerical scheme and the setup of the simulations. We present our main results in § 3, and summarize our conclusions in § 4.

METHOD
To derive the χ eff distribution at mergers, the evolution of the dimensionless BH spins (a) and the binary orbital angular momentum directions (Ĵ bin ) 1 needs to be followed since χ eff is the sum of the projection of massweighted binary spins onto the binary orbital angular momentum, Here m 1 and m 2 are the masses and a 1 and a 2 are the spins of the binary components. To model the evolution of the BH spins and the binary orbital angular momenta, we perform one-dimensional N -body simulations combined with semi-analytic prescriptions. In the following sections, we first give a brief overview of our model and then describe its ingredients in more detail.
2.1. Overview of model In this section, we summarize our model, whose details are presented in Paper I. We consider a system describing a galactic nucleus, consisting of the following five components: (1) a central SMBH, (2) a gaseous accretion disk around the SMBH ("AGN disk"), (3) a spherical stellar cluster, (4) a flattened cluster of BHs, and (5) stars and BHs inside the AGN disk, referred to as the "disk stellar" and "disk BH" components. To follow the time-evolution of the BHs in this system, focusing on their capture by the disk, and the formation, evolution, and disruption of BH binaries in the disk, we run one-dimensional N -body simulations combined with a semi-analytical method. We introduce N -body particles representing either single objects or binaries, and for each particle, we follow its radial position from the central SMBH, as well as its radial velocity, together with the evolution of the binaries' separation. The other two spatial directions are followed only statistically. In this paper, we additionally follow the evolution of BH spins ( § 2.2) and the orbital angular momentum directions of binaries ( § 2.3).
For the AGN disk, we employ the model proposed by Thompson et al. (2005), as adopted in the earlier work by Stone et al. (2017). This represents a Shakura-Sunyaev α-disk with a constant viscosity parameter α and accretion rate in the region where it is not self-gravitating. The model describes a radiatively efficient, geometrically thin, and optically thick disk and extends the disk to pc scales with a constant Toomre parameter in the selfgravitating regime, assuming that it is heated and stabilized by radiation pressure and supernovae from in-situ star formation. We assume that stars and BHs form in the disk at the rate required to stabilize the AGN disk, and some fraction of BHs are initially formed in binaries (see parameter settings in Table 1).
We assume that the AGN disk is surrounded by a spherically symmetric star cluster, with a total mass ≈ 3 times that of the central SMBH within ∼ 3 pc, and with a density profile matching those of a nuclear cluster observed in the Galactic center. We further include a flattened BH cluster component, which has a steeper density profile and a smaller velocity dispersion compared with those of a spherically symmetric star cluster due to mass segregation (e.g. Hopman & Alexander 2006;Szolgyen & Kocsis 2018).
Our model tracks the properties of the BH population, including physical processes due both to the presence of gas and to multi-body dynamical interactions, as follows.
For the interaction with gas, the velocities of all BHs relative to the local AGN disk decrease due to accretion torque and to gas dynamical friction. For binaries of stellar-mass BHs, the binary separation evolves due to gas dynamical friction from the AGN disk and to type I/II migration torque from a small circumbinary disk that forms within the Hill sphere of the binary. Binaries efficiently form in the disk due to gas dynamical friction during two-body encounters (a process we dubbed "gas-capture binary formation"). The radial positions of BHs are also allowed to evolve due to type I/II torques from the AGN disk. Gas accretion affects BH spins and the orbital angular momentum directions of binaries according to newly added prescriptions ( § 2.2.2 and § 2.3.2, respectively).
We also account for dynamical interactions with single stars and BHs and BH binaries. The binaries' separations and velocities evolve due to binary-single interactions, and the velocities of all BHs additionally evolve due to scattering. The evolution of the orbital angular momentum directions of binaries during binary-single interactions are additionally followed with newly added prescriptions ( § 2.3.3). Binaries form due to three-body encounters, and are disrupted by soft binary-single interactions. We also account for GW emission from binaries, which reduces their separation rapidly once they are sufficiently tight. For simplicity, the eccentricity evolution is ignored and orbits around the SMBH and binary orbits are both assumed to be circular.
The interested reader is encouraged to consult Paper I for detailed descriptions of the above model and its ingredients. In the following sections, we only describe the new prescriptions which we added to Paper I, in order to follow the evolution of χ eff of each BH binary.
2.2. BH spin evolution BH spin is characterized by the dimensionless spin parameter a = cJ BH /Gm 2 BH , where G is the gravitational constant, c is the speed of light, m BH is the mass and J BH is the angular momentum of the BH. In this section, we describe the initial distribution and the subsequent evolution of BH spins.

Initial BH spin distribution
In our model, there are two types of BHs differentiated by their origin: BHs formed before the beginning of the current AGN phase (pre-existing BHs), and BHs formed during the current AGN phase (in-situ formed BHs). Pre-existing BHs are distributed in nuclear star clusters, but have a density profile that is steeper (e.g. Hopman & Alexander 2006;Freitag et al. 2006;Keshet et al. 2009) and velocity dispersion that is smaller (Szolgyen & Kocsis 2018) compared to those of typical-mass stars. Pre-existing BHs are expected to have random spin directions since they presumably formed by the disruption of globular clusters (e.g. Mapelli & Gualandris 2016), or by the fragmentation of previous AGN disks or disks in non-active phases whose orientations differed from the current one.
On the other hand, in-situ BHs form in the outer regions of the AGN disk, and could have their spins directed along the angular momentum of the AGN disk. This may be expected if these BHs form from, or efficiently accrete, gas whose angular momentum direction is the same as that of the background AGN disk. This assertion might be justified by the analogy with the planets in the Solar system, all of which except for Venus and Uranus spin in the same direction as their orbital motion to within 30 • .
The typical values of the initial BH spins are highly uncertain. The progenitors of some BBHs, BHs in highmass X-ray binaries, are observed to have high spin (see, e.g. Miller & Miller 2015 for a review). However, we do not have any information on the spins of isolated single BHs or for heavier BHs with masses similar to those discovered by GW observations. We therefore consider several distributions for the initial BH spin a. For the direction of the initial BH spin a, we examine two models: (i) the spin directionâ is random, (ii)â is directed along the angular momentum of the AGN diskĴ AGN (the latter defined with respect to the SMBH), where we fixĴ AGN =ẑ, i.e. along the z-axis. For the magnitude of the initial BH spin |a|, we examine the full range of values between 0 and 0.99 (models M1-M7; see Table 2 below). In the fiducial model (M1), a 0 = 0 for both pre-existing and in-situ BHs. In models M2-M7, we instead assume a 0 = 0.1-0.99, respectively. In all six models M2-M7, the spin direction of the pre-existing and in-situ formed BHs were assumed random and aligned with the AGN disk, respectively. In models M8 and M9, we adopt a 0 = 0.7 for all BHs (e.g. Shibata & Shapiro 2002), but assume that all BH spins are random (model M8), or aligned with the AGN disk (model M9).
Note that in our models the initial BH mass is below 15 M , which may be expected for high-metallicity environments as in AGN. Unlike in other BH formation channels with heavier BHs (e.g. Gerosa et al. 2018), at these lower masses there is no apriori anti-correlation between mass and spin. The BH masses and spins change significantly from their initial values during the evolution in AGN due to gas accretion and mergers in our simulations.

Gas accretion
In our model, BHs capture gas from the AGN disk while they are moving in or crossing the disk. Some fraction of captured gas is assumed to accrete onto BHs through circum-BH disks (for single BHs) or mini-disks (for BBHs, fed from circumbinary disks). During such accretion processes the binary mass, velocity, and separation, as well as the BH spins all evolve.
The spin values a = 1 and −1 represent a maximally spinning BH, and the sign of a is defined so that for a > 0 the BH is spinning in the same direction as the inner accretion disk, and for a < 0 the spin is in the opposite direction. The spin magnitude after an accretion episode is given by (2) (Bardeen 1970), where f acc ≡ (m BH + ∆m BH )/m BH , ∆m BH ≡ṁ BH ∆t is the mass accreted during the time step ∆t, the superscript f stands for the values after the episode, and r isco is the radius of the innermost stable circular orbit (ISCO) in reduced units, defined as (3) with the minus sign for a > 0 and the plus sign for a < 0. Here R g = Gm BH /c 2 is the gravitational radius of the BH, and the functions Z 1 and Z 2 are given by While Eq.
(2) gives unphysical spin values when a highly spinning BH accretes a significant amount of gas (a f > 1 or imaginary), the torque exerted by the radiation of a thin accretion disk prevents a f > 0.998 (Thorne 1974). Fully relativistic magnetohydrodynamics simulations suggest that the spin value does not grow beyond 0.95 during accretion from a thick disk (Gammie et al. 2004;Shapiro 2005). We set the upper limit of a f to 0.99 in our models. When the spin angular momentum J BH = a Gm 3 BH R g of a BH is misaligned with its inner disk, the BH induces a Lense-Thirring precession in the misaligned disc elements, which causes the inner parts of the disk and the BH spin to align. The transition between aligned and misaligned annuli of the disk occurs at the so-called warp radius R warp . In each time-step ∆t in our model, J BH aligns with the initial spin angular momentum of the BH plus the angular momentum ∆J warp of the disk within the warp radius: J BH → J BH + ∆J warp . For a Shakura-Sunyaev disk, the warp radius is given by (e.g. Volonteri et al. 2007), where f Edd =ṁ BH c 2 /L Edd is the accretion rate in Eddington units (without a radiative efficiency), L Edd is the Eddington luminosity, ν 1 is the viscosity responsible for transferring angular momentum in the accretion disk, and ν 2 is the viscosity responsible for warp propagation. We set ∆J warp = ∆M warp Gm BH R warp (e.g. Volonteri et al. 2007), where ∆M warp is the warped disk mass aligned with the BH spin during ∆t. We set ∆M warp = ∆m BH , and ∆Ĵ warp points in the same direction as the angular momentum of the circum-BH diskĴ CBHD . After alignment, the magnitude of the BH spin evolves through gas accretion via Eq.
(2). We note that whenever the condition is satisfied, where θ BH,warp is the angle between J BH and ∆J warp , the accretion disk within R warp becomes antialigned with the BH spin direction (King et al. 2005). Due to the Lense-Thirring effect, the BH spin and the circum-BH disk angular momentumĴ CBHD align faster than how the magnitude of the BH spin grows (e.g. Volonteri et al. 2007). In this process, there are two large uncertainties: the size of the warp radius, and the direction ofĴ CBHD The size of the warp radius strongly depends on the ratio ν 2 /ν 1 (Eq. 6). Many studies adopted ν 2 /ν 1 = 2(1 + 7α SS )/(4 + α 2 SS )/α 2 SS ∼ 85, motivated by analyses of low-amplitude warps (Ogilvie 1999). On the other hand, Lodato & Gerosa (2013) have shown that when considering large-amplitude warps, ν 2 /ν 1 is between ∼ 2 − 50 depending on α SS and the misalignment angle θ BH,warp . We set ν 2 /ν 1 to be a free parameter with a fiducial value of 10, and vary it from 2 to 50 (models M10 and M11).
For a single BH,Ĵ CBHD aligns withĴ AGN (see results in Lubow et al. 1999 in the context of protoplanetary disks). On the other hand, when a BH is in a binary, J CBHD aligns with the orbital angular momentum direction of the binaryĴ bin on the disk's viscous timescale (e.g. Moody et al. 2019). Assuming a Shakura-Sunyaev disk, the viscous timescale is given by (e.g. Frank et al. 2002), where s is the binary separation, m bin is the binary mass,ṁ bin is the accretion rate onto the binary,ṁ Edd,bin = L Edd /(η c c 2 ) is the Eddington accretion rate for the binary, and η c is the radiative efficiency. Due to the short viscous timescale, in our fiducial setting we assume thatĴ CBHD is the same asĴ AGN when a BH is single, and is the same asĴ bin when a BH is in a binary. For comparison, we also investigate the alternative assumption that the direction ofĴ CBHD is always aligned withĴ AGN (model M12). We set the BH accretion rate to the minimum of the Eddington accretion rate and the Bondi-Hoyle-Lyttleton rate (Eq. 24 in Paper I) taking into account a reduction due to the shearing motion of the nearby disk gas (Eqs. 29-32 in Paper I). When BHs are in binaries, we apportion the total accretion between primary and secondary BHs following Duffell et al. (2019), which is updated from Paper I in which we used earlier results from Farris et al. (2014).

Mergers
Following BBH mergers, the dimensionless spin parameter a f of the remnant BH depends on the spins a 1 and a 2 and the dimensionless orbital angular momentum parameter l of the two original binary components. 2 We adopt the formula obtained from numerical simulations of BBH mergers in Rezzolla et al. (2008), where q ≡ m 2 /m 1 ≤ 1 is the mass ratio. The magnitude of l is given by where η ≡ q/(q + 1) 2 is the symmetric mass ratio, s 4 = −0.129, s 5 = −0.384, t 0 = −2.686, t 2 = −3.454, and t 3 = 2.353 are values obtained in Rezzolla et al. (2008), and θ 12 , θ 1b , and θ 2b are the angles between the spins of the two BHs and their orbital angular momentum, According to Eq. (9), if a 1 = a 2 = 0, a f monotonically increases from 0.58 to 0.69 as q increases from 0.4 to 1 and decreases to 0 as q → 0. Following Rezzolla et al. (2008) and Dubois et al. (2014), we assume that GW radiation does not affect the direction of the orbital angular momentum, and set the direction of l to the binary orbital angular momentum at merger. Although Barausse & Rezzolla (2009) suggested that GW radiation modifies the binary orbital angular momentum direction just before merger, we neglect this correction for simplicity as results are not sensitive to the precise a f direction. Similarly, we assume that χ eff values are unaffected by relativitistic effects. This is justified if the BH spin directions and the orbital angular momentum directions are not influenced by the effects before binaries enter the frequency above which the LIGO/Virgo detectors are sensitive, which is ∼ 10 Hz. The orbital angular momentum direction is not directly influenced by GW radiation reaction significantly for circular orbits at lower frequencies (Barausse & Rezzolla 2009). The BH spin directions can be systematically affected by spin-orbit resonances due to precession of the BH spins, which can align or anti-align the BH spins with each other or cause nutation (Kesden et al. 2010;Gerosa et al. 2019). In our simulations, such resonances or nutation occur only within the detectable frequency band above 10 Hz for 95% of stellar BH mergers. Thus, we conclude that these general relativistic effects have a small impact on the detectable χ eff distribution for LIGO/VIRGO. 2 (G/c)m 1 m 2 l is the orbital angular momentum at the ISCO 2.3. Evolution of binary orbital angular momentum direction 2.3.1. Initial orbital angular momentum direction We consider four types of BH binaries distinguished by their formation process: (i) pre-existing binaries, (ii) gascapture binaries, (iii) dynamically formed binaries, and (iv) remnants of stellar binaries formed in-situ. For binaries belonging to (iii) we setĴ bin = ±Ĵ AGN , and the ratio of aligned binaries over anti-aligned binaries to 1, as suggested by simulations of binary formation in migration traps (Secunda et al. in prep). For binaries belonging to (iv), we also set the ratio to be 1 for simplicity, although we note that this ratio is highly uncertain. On the other hand, for binaries belonging to (i) or (ii), the orbital angular momentum directions are presumed to be random. For simplicity, we assume that all binaries have zero eccentricity. We expect that this assumption does not significantly change the evolution of binary orbital angular momenta or BH spins.
As we show below, the initial angular momentum directions of binaries have a relatively small effect on the χ eff distribution measured at merger, because these directions are frequently randomized by binary-single interactions.

Gas accretion
The orbital angular momentum directions of binaries evolve due to accretion torques. If the circumbinary gas is rotating in the same direction as the AGN disk, the binaries are aligned with the AGN disk since the relative velocity between the binary components and the gas is reduced by the accretion. Referring to Lubow et al. (1999), we assume that the angular momentum direction of the captured gas with respect to the binary (Ĵ gas ) is the same as the angular momentum direction of the AGN disk with respect to the SMBH (Ĵ AGN ) for binaries embedded in the AGN disk. The angular momentum of the captured gas is added to the binary orbital angular momentum as J f bin = J bin + J gas , where we set and where v bin (s) = Gm bin /s is the relative rotation velocity of binary components, and f rot is a parameter determining the efficiency of the alignment of the binary angular momentum direction due to gas capture. In the fiducial model, we set f rot = 1 assuming that the binary receives a torque from gas circularly rotating at ∼ s from the binary. However, a low degree of rotation (f rot ∼ 0) of gas accreting onto a low-mass object in an AGN disk is suggested in simulations by Baruteau et al. (2011) and Derdzinski et al. (2018). For completeness, we investigate this case, as well as an opposite extreme case with f rot = 10 (models M13 and M14). Note that during gas accretion, the binary separation also evolves due to type I/II torque of the circumbinary disk (Paper I).
2.3.3. Binary-single interaction After a hard binary-single interaction, the orbital angular momenta of binaries are modified due to chaotic interactions. In this paper, we simply assume that after a hard binary-single interaction, the orbital angular momentum direction of a binary becomes isotropically random.
Whenever a binary-single interaction occurs in the simulation, we choose a nearby third object, and assign a recoil kick velocity to it. If this third object is itself a binary, we assume that the softer binary is disrupted, while the harder binary experiences a hard binary-single interaction by regarding the softer binary as a single object for simplicity, and assign the recoil kick velocity to its center of mass.

Merger prescription
Since we track the evolution of the BH spins and the binary orbital angular momenta, we can estimate the recoil velocity due to anisotropic GW radiation and mass loss at mergers more precisely. We add the following prescriptions to the model used in Paper I.

Recoil velocity at merger
Due to the burst of anisotropic GW radiation at merger, a remnant BH receives a recoil kick. To calculate the recoil velocities, we adopt the fitting formulae obtained from numerical simulations by Lousto et al. (2012), where v m is a mass-asymmetry contribution, v ⊥ and v are kick components perpendicular and parallel to the orbital angular momentum, respectively, e x ,ê y are orthogonal unit vectors in the orbital plane, andê z is the direction of the binary orbital angular momentum. The symbols and ⊥ refer to the directions parallel and perpendicular to the orbital angular momentum, respectively, and the numerical constants are A = 1.2 × 10 4 km/s, B = −0.93, H = 6.9 × 10 3 km/s, ξ = 145 • ± 5 • , V 1,1 = 3678 km/s, V A = 2481 km/s, V B = 1792 km/s, and V C = 1507 km/s, φ 1 is the phase angle of the binary, and φ ∆ is the angle between the in-plane component of the vector and the infall direction at merger. We choose φ ∆ − φ 1 from a random distribution uniform in [0, π].
2.4.2. Mass loss at merger Due to GW radiation, some fraction of the BH mass is radiated away. We adopt a simplified approximation for the remnant mass from Barausse et al. (2012), whereẼ is the energy per unit mass of a particle with spiñ and p 0 = 0.04827 and p 1 = 0.01707 are parameters obtained by fitting numerical results.

Merger condition
We assume that a binary merges when its separation s becomes smaller than the ISCO of a particle with mass m bin and spin given in Eq. (20). We subsequently treat the object as a single BH with a mass given by Eq. (18). Table 1 lists the parameter values adopted in the fiducial model (the same as Model 1 in Paper I). We assume that stars are distributed spherically with a Maxwell-Boltzmann velocity distribution, the stellar mass within 3 pc is M star,3pc = 10 7 M , and the power-law slope of the stellar initial mass function (IMF) is δ IMF = 2.35. BHs are initially distributed from r in,BH = 2 × 10 −4 pc to r out,BH = 3 pc with a cumulative radial profile

Numerical choices
with a power-law index γ ρ = 0, where N BH,ini (r) labels the total initial number of BHs within a distance r from the central SMBH. Note that BHs are assumed to have a flattened axisymmetric distribution (see Paper I). The x, y, and z velocities for BHs relative to the local Keplerian value v Kep (r) are initially drawn from a Gaussian distribution with the dispersion of β v v Kep (r)/ √ 3, where β v is a velocity anisotropy parameter set to β v = 0.2 motivated by vector resonant relaxation (Szolgyen & Kocsis 2018). The total number and mass in BHs are calculated from the stellar mass, the stellar IMF, and the relation between the stellar and BH mass derived in Belczynski et al. (2010). We set the fraction of pre-existing binaries to be f pre = 0.15. In the fiducial model, there are initially 2 × 10 4 BHs and 1.5 × 10 3 binaries. As in Paper I, the time step parameter is η t = 0.1, and the number of radial cells storing physical quantities is N cell = 120.
The mass of the central SMBH is M SMBH = 4×10 6 M , the gas accretion rate from the outer radius isṀ out = 0.1Ṁ Edd , whereṀ Edd = L Edd /(η c c 2 ) is the Eddington accretion rate, here defined including a radiative efficiency of η c = 0.1. The efficiency of angular momentum transport in the α-disk is α SS = 0.1 (Shakura & Sunyaev 1973), and the angular momentum transfer parameter is m AM = 0.15 (Thompson et al. 2005).
3. RESULTS 3.1. χ eff evolution: an illustrative example We used the combination of semi-analytical calculations and simulations, described above, to investigate the χ eff distribution of BHs merging in AGN disks. An illustrative example of the evolution of such a BH binary is shown in Figure 1. Table 1 Fiducial values of our model parameters.

Parameter
Fiducial value Initial BH spin magnitude |a| = 0 Angular momentum directions of circum-BH disksĴ CBHD =Ĵ AGN for single BHs, J CBHD =Ĵ bin for BHs in binaries Ratio of viscous parameters ν 2 /ν 1 = 10 Efficiency of alignment of J bin due to gas capture frot = 1 Mass of the central SMBH M SMBH = 4 × 10 6 M Gas accretion rate at the outer radiusṀout = 0.1Ṁ Edd Fraction of pre-existing binaries fpre = 0.15 Power-law exponent for the initial density profile for BHs γρ = 0 Parameter setting the initial velocity anisotropy for BHs β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 disk m AM = 0.15 Accretion rate in Eddington units onto stellar-mass BHs Γ 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 The binary in this figure forms via the gas-capture mechanism at t = 3.5 Myr at a distance of 1.3 pc from the SMBH, with an initial separation of 3.4 × 10 −3 pc. The masses of the binary components are 10.7 and 6.8 M , the magnitude of BH spins are 0.79 and 0.70, and the angles between the BH spins andĴ AGN are 0.20 and 3.1, respectively. 3 The separation of the binary (black line in panel (a)) decreases, successively, due to gas dynamical friction, binary-single interactions, and GW radiation as shown in Paper I.
While the binary is in the AGN disk (v z < c s , blue and orange lines in panel (b)),Ĵ bin aligns withĴ AGN due to accretion torque (orange line in panel (d)). Also,â 1 and a 2 (blue and cyan lines in panel (d)) evolve towards the angular momentum direction of a circum-BH disk, which is set to be the same asĴ bin . Such alignment ofâ 1 andâ 2 withĴ bin increases χ eff (panel (c)). The spin magnitudes, |a 1 | and |a 2 | evolve due to gas accretion (black and gray lines in panel (d)), but only by 20% and 11%, which are much smaller than the change in χ eff . Until t ∼ 5 Myr, since the anti-alignment condition (Eq. 7) is satisfied for the secondary BH, |a 2 | slightly decreases as gas accretes.
After each binary-single interaction,Ĵ bin is randomized, which reduces |χ eff | on average. This binary merges outside the AGN disk 17.6 Myr after it formed, and its components accrete 2.0 and 1.7 M until their merger. Sinceâ 1 andâ 2 are roughly anti-aligned withĴ bin following binary-single interactions, χ eff at merger has a negative value of -0.59. Thus, χ eff in this case evolves through both gas accretion and binary-single interactions.

χ eff distribution in different models
In this section, we present the probability distribution of χ eff as well as of several different quantities at the time of merger, obtained in our models. Because of the flexibility of our model, we are also able to study the dependence of these results on the choice of prescriptions and parameter values. In Table 2, we list the model variations we have investigated. These include the fiducial model (M1), and 23 different variations (models M2-M24). The variations can be divided into four categories. First, we examine different choices of the initial BH spin magnitudes and directions (models M2-M9). Second, we vary the prescriptions for the evolution of the BH spin directions and the binary angular momentum directions during gas accretion (models M10-M14). Third, we study different parameters related to the AGN disk (models M15-M19). Finally, we vary the properties of the initial BH population (models M20-M24). Fig. 2 shows our main results, namely the differential probability distributions of χ eff , |a|, cos θ a , cos θ bin , and m chirp (first to fifth column, respectively). Here, |a| stands for either |a 1 | or |a 2 | and θ a and θ bin are the angles between (â) and the binary orbit (Ĵ bin ) with respect to the AGN disk (Ĵ AGN ). The four different rows in this figure correspond to the four different types of model variations, as discussed above, and labeled in the figure.
The distributions of all predicted quantities in Fig. 2 are weighted by the detectable volume, which enables us to compare our predictions to the observed distribution (dashed red lines). The detectable volume is calculated via Eq. (6) of The LIGO Scientific Collaboration & The Virgo Collaboration (2012) and using the noise spectral density of the ER13 (prior to O3) run of LIGO Hanford (Kissel & Betzwieser 2018), in which the volume is assumed to depend only on the masses of the binary components. The volume is roughly proportional to m 2.2 1 for m bin 100 M (e.g. Fishbach & Holz 2017). Furthermore, to compare with the observed distributions of χ eff and m chirp , we add observational errors to χ eff and m chirp for each merger in the simulations. 4 For a simple treatment, we draw the errors of χ eff and m chirp from independent Gaussian distributions whose 90% intervals are ±0.2 and ±0.08m chirp , respectively, which match the typical O1/O2 observational error magnitudes (The LIGO Scientific Collaboration et al. 2018). These errors are added to our predictions for the analysis of the Kolmogorov-Smirnov (KS) test and Figs. 2-9 below.
In Fig. 3 we present the cumulative, rather than the differential probability distributions shown in Fig. 2. We present and discuss results at t = 10 Myr unless stated otherwise (but see Table 2 below for different choices).
The fiducial model is shown by black lines in all of the panels, (a)-(e), of Figs. 2 and 3. Panel (d) shows that θ bin represents an isotropic distribution (uniform in cos θ bin ). This is becauseĴ bin is frequently randomized by binarysingle interactions. On the other hand,â 1 andâ 2 tend to align withĴ AGN (small θ a ) for the following reason. First,â 1 andâ 2 are gradually aligned withĴ bin due to gas accretion. HoweverĴ bin aligns withĴ AGN when the binary is in the AGN disk and it is mostly random when outside of it. Thus, the θ a distribution is influenced by the fraction of the time that binaries typically spend in the AGN disk.
The mean spin magnitude |a| evolves from 0 in the fiducial model to typically lie in the range 0.08-0.64, χ eff and m chirp are typically in the range −0.22 to +0.24 (−0.18 to +0.21 if no errors are added) and 8-49 M , respectively (enclosing 68% of the total probability). By comparison, for first-generation 5 mergers |a|, χ eff , and m chirp are distributed in the range 0.044-0.32, −0.15 to +0.20, and 7-11 M , respectively. Hence, |a| and m chirp both evolve significantly due to mergers. This trend is also visible in Fig. 11, which shows that the standard deviation of χ eff increases with m chirp (orange line in panel (a)). It is notable that a similar trend is also seen in the observed distribution (orange line in panel (t)).
To see the robustness of this trend, in Fig. 5 we plot the weighted mean and standard deviation of χ eff as a function of mass for ten additional simulations for independent realizations of the initial condition. The orange line in Fig. 5 Table 2. quantity x i for each model is calculated by where i is the index of a merger, w i is the detectable volume, j is the index of a bin, M j is the number of nonzero weights, andx j = i∈j w i x i / i∈j w i is the weighted mean of x i in the j th bin. In Fig. 5, the errors of the mean and the standard deviation are calculated by and σ(s j ) = s j 1 2(M j − 1) respectively (Harding et al. 2014). Eq. (24) is approximately correct for M j 10. Fig. 5 shows that the standard deviation of χ eff robustly increases up to m chirp ∼ 20 M , while it is roughly constant in the range 20 m chirp 100.
Next, we present the dependence of the distributions on the assumed initial BH spins (first row in Fig. 2). In models M1 and M7, the initial BH spins are set to |a 0 | = 0 and 0.99, respectively. We also examined five intermediate cases, with |a 0 | = 0.1, 0.2, 0.3, 0.5, and 0.7 (models M2-M6) and found that the resulting distributions for all five quantities lie in-between those of the extreme models M1 and M7. For clarity, the intermediate cases (M2-M5) are therefore not shown in Figs. 2 and 3. Gas accretion typically does not cause a systematic shift, but rather smears the initial distribution of |a| (dashed line in panel (b) of Fig. 3). It typically increases the spin along the orbital angular momentum which is frequently reoriented in random directions by binary-single interactions). In contrast, |a| evolves to ∼ 0.7 after mergers (solid and dashed black and blue lines). Models M1 and M7 show that |a 0 | has an influence on the |a| and the χ eff distribution. This is seen as the difference between the black and the cyan lines in panels (a) and (b) in Figs. 2 and 3. This suggests that the initial spin distribution might be constrained by GW observations if they originate in AGN disks.
To examine the effects of the initial spin directions, in model M8 we set |a 0 | = 0.7 and draw randomâ 0 directions and in model M9 we set a 0 = 0.7Ĵ AGN for all (i.e., including pre-existing) BHs. The resulting distributions for model M8 are similar to those in model M6 (orange and brown lines in the top row). This is because the difference between models M6 and M8 is only the initial spin directions for in-situ formed BHs, whose contribution to all mergers is small (∼ 1%). The alignment of a 0 withĴ AGN in model M9 might be realized if previous AGN episodes yield spin directions aligned with the present-day disk orientation. We include this extreme model as an academic exercise to investigate the impact of full initial spin alignment. In this model, M9, θ a is distributed around low values due to the initial direction of a 0 , while M6 and M8 have broader distributions. The orange and green lines in panel (a) suggest that the a 0 direction has a negligible influence on the χ eff distribution. Models M10-M14 (second row in Figs. 2 and 3; panels (f)-(j)) show the results when changing the prescriptions or the parameters affecting the evolution ofâ and J bin during the accretion episodes. When the BH spins align efficiently due to enhanced viscosity accelerating the Bardeen-Petterson effect (model M11) or when the angular momentum of gas captured by the binaries is random (model M13), θ a is slightly closer to an isotropic distribution (orange and brown lines in panel (h)). Nevertheless, we can see that these changes have a small impact on the χ eff distribution (panel (f)). This is simply because the timescale of randomization ofĴ bin due to binary-single interactions is shorter than the timescale of alignment ofâ towardĴ bin due to gas accretion. Likewise, these prescriptions have very little impact on the distribution of the other quantities.
In the third rows of Figs. 2 and 3, we examine the impact of AGN disk properties. It is not clear whether (or to what extent) radial migration operates due to the complexity of the effects of N -body migrators (Broz et al. 2018), feedback from BHs (del Valle & Volonteri 2018; Regan et al. 2019), and inhomogeneities in the turbulent accretion disk (Laughlin et al. 2004;Baruteau & Lin 2010). In model M15 (shown in black lines), radial migration due to torques from the AGN disk is assumed to be inefficient. In this model, θ a and θ bin (panels (m) and (n)) are distributed around lower values compared with those in the fiducial model. This is because binaries cannot migrate to high BH-density regions, where disorienting binary-single interactions are frequent (Paper I). As a result, χ eff is distributed toward higher values (panel (k)). Also, m chirp tends to be lower since repeated "hierarchical" mergers, which build up the more massive BHs, are less frequent (panel (o)). We note here that except for model M15, the m chirp distributions are very similar in all model variants (panels (e), (j), (o), and (t)). This is because the m chirp distributions are determined primarily by how often repeated mergers occur,  The probability that the χ eff and m chirp distributions for all events is reproduced by each model (the AGN contribution to all merger is f AGN = 1) estimated by the KS test (P KS,χ eff ,m chirp ) is shown in a lower right corner. Panel (t) shows the distribution derived by performing a kernel-density estimate for the observed distribution, which is presented by different colors to emphasize its peculiarity. The values for the events inferred from the LIGO/Virgo collaboration (red circles) and the IAS group (cyan circles) are overplotted in all panels. and this is not significantly influenced by the parameters we changed, other than the efficiency of migration. On the other hand, the χ eff and m chirp distributions for a model without migration would be sensitive to various other parameters, such as the initial mass distribution of BHs, and the AGN lifetime.
In the bottom row of Figs. 2 and 3, in models M16-M24, we examine the influence of the parameters of the AGN disk or the initial BH distribution (see Table 2). The resulting χ eff distributions are similar in these models (panels (k) and (p)). This is again because the χ eff distribution is mainly affected by how frequently spindisorienting binary-single interactions take place, which is not sensitive to the changes mentioned above.
Overall, we find that the χ eff distribution is relatively sensitive to the values of |a 0 | and the efficiency of migration, and the m chirp distribution is significantly influenced by migration. In the next section, we compare these predictions with observations.

Comparison with observed distribution
In Figs. 6 and 7, we compare the χ eff distributions predicted by our models with that inferred from the observed GW events. In the observed distribution reported by the LIGO and Virgo collaborations (The LIGO Scientific Collaboration et al. 2018) the χ eff values are concentrated at low absolute values, near zero (Fig. 7). On the other hand, a few possible additional events have been identified with higher and lower χ eff values (Fig. 6, Zackay et al. 2019a,b;Venumadhav et al. 2019; but see also Huang et al. 2020). Note that the IAS group (e.g. Zackay et al. 2019b) and the Hannover group (Nitz et al. 2020) also recovered the events reported by the LIGO/Virgo collaborations.
In Fig. 4, we additionally show the detection rate distributions predicted in several models in the χ eff vs. m chirp plane, together with the distribution inferred from the observed GW events (Table 3).
To compare the predicted and the observed distributions quantitatively, we use the KS test as well as a Comparison between the χ eff distribution inferred from the observed GW events and those predicted in our models. Black, orange, cyan, and brown lines show the distribution in models M1, M4, M7, and M15, respectively. The predicted distributions are weighted by the detectable volume, and include observational errors. Red and blue circles show the median χ eff values reported by The LIGO Scientific Collaboration et al. (2018) and the IAS group (Zackay et al. 2019b,a;Venumadhav et al. 2019), respectively. Error bars correspond to 90% credible intervals.
Bayesian analysis. The KS test enables us to estimate the probability that the distribution of all or a subset of the observed events is reproduced by a given model, while the Bayesian analysis can be used to assess how consistent each individual event is with a given model. Table 2 lists the results of the KS test. In each model, P KS,χ eff , P KS,m chirp , and P KS,χ eff ,m chirp are the probabilities that the set of all measured χ eff and m chirp values were drawn from the one-dimensional χ eff , m chirp , and the joint two-dimensional (χ eff , m chirp )-distributions  Figure 8. The KS probabilities that the one-dimensional χ eff distribution inferred from observed GW events is consistent with models assuming different initial spin magnitudes |a 0 | (i.e. models M1-M7). Black and orange lines show the results in which all observed GW events and only the LIGO/Virgo events are used, respectively. Solid and dashed lines show the results in which errors are and are not included in the predicted χ eff distribution, respectively. For a 0 = 0, the errors and means are calculated by performing ten additional runs with different realizations of the initial conditions. The observed χ eff distribution slightly favors moderate values of |a 0 | 0.5. predicted in that model, computed following Press & Teukolsky (1988). P KS,LV,χ eff and P KS,LV,m chirp are the probabilities that χ eff and m chirp values for the events reported by the LIGO/Virgo collaborations were drawn from the predicted one-dimensional χ eff and m chirpdistributions, respectively. Fig. 4 also lists P KS,χ eff ,m chirp in each panel. As described in Figs. 2 and 3, each predicted merger is weighted by the detectable volume, and errors are added on the predicted χ eff and m chirp values. The values of P KS,χ eff ,m chirp are typically ∼ 0.01-0.1, except for model M15, which yields a much lower value of P KS,χ eff ,m chirp = 6.5 × 10 −6 . This is because m chirp is typically much lower in model M15 compared to the other models, as well as compared to the observations (panel (k) in Fig. 4, P KS,m chirp = 9.2 × 10 −9 ). As explained above, this is because in this model (M15), radial migration is turned off; this makes hierarchical mergers much less common. We note, however, that the m chirp distribution has large uncertainties in our models, for several reasons. First, we do not take into account the exchange of binary components during binary-single interactions, which affects the m chirp distribution (the main assumptions in our models are listed in § 2 of Paper I). Also, the time evolution of the AGN disk model, which we ignore, may affect the merged mass distribution ( § 5.7.1 of Paper I). Indeed, using a 30 Myr AGN lifetime in Paper I we found that the mass distribution of merging BHs extends to the values matching the observations if radial migration is turned off (see Figure 14 therein, panels a and d). Additional uncertainties include the stellar IMF in galactic centers (Lu et al. 2013) and the relation between the initial stellar mass and its remnant BH mass (e.g. Belczynski et al. 2010;Chen et al. 2015). Particularly, we neglected the possibility that BHs may be delivered to the nuclear star cluster by the infall of low metallicity globular clusters where the BH masses are expected to be higher (Tremaine et al. 1975;Antonini 2013;Arca-Sedda et al. 2018; Arca Sedda & Benacquista 2019; Arca Sedda 2020). For a more rigorous comparison with the observed m chirp distribution, these points should be considered in a future study.

Kolmogorov-Smirnov test
Focusing only on the χ eff distribution, we find P KS,χ eff are ∼ 0.1 − 0.7. This suggests that the observed χ eff distribution is consistent with most of our models. Fig. 8 shows P KS,χ eff (black lines) and P KS,LV,χ eff (orange lines) as a function of |a 0 | (models M1-M7). The solid and dashed lines show the results in which observational errors are and are not included, respectively. P KS,χ eff and P KS,LV,χ eff are highest (0.68 and 0.15) in model M3 and M4 in which |a 0 | = 0.2 and |a 0 | = 0.3, respectively. Thus, moderate values for |a 0 | are preferred by the observed χ eff distribution.
We further investigate how the KS probabilities change if we assume that only some fraction f AGN < 1 of mergers occur in AGN disks, with the remaining fraction (1 − f AGN ) produced in other unrelated channel(s). This gives an estimate for the maximum allowed fraction of events related to AGN disks in each of our models. For non-AGN mergers, we conservatively assume that the χ eff and m chirp distributions are the same as the observed distributions to date (Table 3), but the detection rates are normalized to (1 − f AGN ). On the other hand, for AGN mergers, each merger is weighted by the detectable volume referring to Kissel & Betzwieser (2018) as before, and the total detection rate distribution of χ eff and m chirp in each model is normalized to f AGN .
We construct the χ eff and/or m chirp distributions by summing non-AGN and AGN mergers, and calculating the KS probability of their combined distribution. The cyan, orange, and black lines in Fig. 9 show P KS,χ eff , P KS,m chirp , and P KS,χ eff ,m chirp as a function of f AGN . The thick solid lines are the probabilities including all events (Table 3), and the dashed lines include only the events reported by the LIGO/Virgo groups (The LIGO Scientific Collaboration et al. 2018). Even for model M15, in which P KS,χ eff ,m chirp is lowest for f AGN = 1, we find P KS,χ eff ,m chirp is 0.7 provided that f AGN 0.30 (solid black line in the right panel). Thus, at least ∼ 30% of mergers might originate in AGN disks even in the worst model (see discussion above on caveats which may increase f AGN for model M15).

Bayesian analysis
Next, to assess the relative likelihood to produce each event in different models, we calculate the Bayes factors between pairs of models, where P (d i |A) is the likelihood of obtaining a data d i in an event i from Model A, where P (d i |m chirp , χ eff , q) is the three dimensional likelihood for m chirp , χ eff , and q, and P (m chirp , χ eff , q|A) is the probability distribution of m chirp , χ eff , and q in Model A.
To calculate P (m chirp , χ eff , q|A), we first count mergers in 30 × 30 × 30 uniform bins in χ eff , m chirp , and q for Model A. The maximum and minimum values of m chirp for the bins are set to 150 and 5 M , respectively. In this procedure, we weighted each merger by the detectable volume. To reduce the statistical fluctuation in the distribution of χ eff , m chirp , and q due to the finite number of mergers in our models, we perform a kernel-density estimate for the distribution using Gaussian kernels whose bandwidth is chosen to satisfy Scott's Rule (Scott 1992).
For simplicity, we assume that they follow independent Gaussian distributions, as commonly assumed in studies analyzing observed GW data (e.g. Fishbach & Holz 2017): where N (c 1 , c 2 2 ) is the Gaussian distribution with average c 1 and dispersion c 2 2 , m chirp,i , χ eff,i and q i are the median values and σ 2 m chirp ,i , σ 2 χ eff ,i , and σ 2 q,i are the dispersions of m chirp , χ eff , and q observed in GW event i. We set the average values and the standard deviation for each event to the median values and the 90% credible intervals in Table 3 divided by 3.3, as appropriate for a Gaussian distribution. For the events found by the IAS group, for simplicity, we calculate the dispersion of the source-frame chirp mass assuming no covariance between the parameters.
We calculate P (d i |m chirp , χ eff , q) by generating 1000 samples according to the Gaussian distribution and normalizing the distribution as The values for P (d i |m chirp , χ eff , q) are stored in 30 × 30 × 30 uniform bins. For each event i, we calculate the Bayes factor for a Model A relative to the observed distribution (i.e. "Model B" in the ratio in Eq. 25 is taken to be the observed distribution itself). The observed distribution is constructed by smoothing the observed m chirp,i , χ eff,i and q i distribution using a kernel density estimate as applied above (panel (t) of Fig. 4). This K A,obs,i presents the strongest test of Model A for each event, since models are compared to the actual observed distribution.
In most models, the lowest value of K A,obs,i among the GW sources is typically ∼ 0.02-0.2 (min i K A,obs,i in Table 2), except for M15, in which it is much lower (∼ 10 −5 ). In model M15, K A,obs,i is lowest for GW170817A, which is the source with the highest m chirp . These findings are consistent with Fig. 4, which shows that the observed value for the event is outside the predicted range. For the fiducial model M1, we list the value of K A,obs,i for each observed source in Table 3. The most constraining event is GW151216 (K A,obs,i = 0.050); the source which has the highest χ eff (see also Huang et al. 2020). Thus, as expected, the events with the highest χ eff and m chirp constrain the models most strongly.

Comparison to other formation channels
In this section, we briefly discuss differences in the expected distributions of χ eff and/or m chirp between the AGN disk-assisted channel and other binary merger channels.
First, our models with low |a 0 | produce the positive correlation between m chirp and the dispersion of χ eff in  Fig. 4). The field binary evolution channels likely favor rather negative correlation between m chirp and the dispersion of χ eff (Gerosa et al. 2018;Bavera et al. 2019;Safarzadeh et al. 2020). The positive correlation is expected for repeated mergers, which frequently occur for multi-body systems in high escapevelocity environments such as galactic nuclei (Arca , and/or if initial BH spins are low (O'Leary et al. 2016). On the other hand, Arca  show that the positive correlation is not reproduced by mergers in dynamical environments. Hence, the positive correlation suggested by Safarzadeh et al. (2020) may be a signature that the observed mergers are facilitated in AGN disks.
Second, AGN disks can produce high-m chirp mergers. For field binaries, m chirp is limited to 40 M due to pair instability supernovae (e.g. Kinugawa et al. 2014;Spera et al. 2019). In the scenarios involving dynamical formation and evolution, Arca  predicted that 99% of mergers have m chirp 50 M . For mergers in AGN disks, we find that ∼ 10-15% of mergers have m chirp 50 M if BHs migrate efficiently. Thus, if mergers with m chirp 50 M are discovered, they would favor the AGN-disk origin. Although a false alarm rate is high (0.34 yr −1 ), Udall et al. (2019) reported a highmass binary BH merger event, GW170502, with a chirp mass of ∼ 70 M . Similar events with high S/N ratio will be an additional signature for mergers in AGN disks.
Third, if migration is inefficient (model M15), χ eff can be negative and the χ eff distribution may lack symmetry around χ eff = 0. From Table 2, the absolute value for the 90 percentiles for χ eff (|χ eff,90 |) is larger than that for the 10 percentiles (|χ eff,10 |) by ∼ 0.27 in this model. Such asymmetric distribution of χ eff is caused by gas accretion, while it is reduced by the randomization of the binary angular momentum directions due to binarysingle interactions.
Mergers in isolated environments are unlikely to produce negative χ eff (Bavera et al. 2019), unless angular momentum transfer is inefficient; however in this case high-χ eff mergers are overproduced . Mergers in globular clusters and galactic nuclei (without AGN disks) can produce negative χ eff , but the χ eff distribution is almost perfectly symmetric around χ eff = 0 (Rodriguez et al. 2018). Hence, if the asymmetric distribution is observed, a possible interpretation is that mergers originate in AGN disks and binary-single interactions are less efficient (e.g. due to inefficient inward migration to the densely populated inner regions).
In summary, the positive correlation between m bin and the dispersion of χ eff , high-m chirp mergers, and an asymmetric χ eff distribution might be possible signatures that distinguish mergers in AGN disks from other channels. However, the χ eff and m chirp distributions for mergers in AGN disks are found to be strongly affected by radial migration of BHs. The efficiency of this migration is still poorly understood, and should be investigated in the future.

Consistency with GW190412
Recently, a low mass-ratio event, GW190412, has been reported (The LIGO Scientific Collaboration & the Virgo Collaboration 2020). This is the first event which has a low mass ratio (q = 0.28 +0.13 −0.07 ), is constrained to have non-zero BH spin parallel to the binary's orbital plane (χ p = 0.30 +0.19 −0.15 ), and has a primary BH with large spin (a 1 = 0.43 +0.16 −0.26 ). Fishbach & Holz (2020) predicted that 99% of mergers have q > 0.51 from the LIGO/Virgo events in O1/O2, which suggests that GW190412 is a highly unusual event. Here, we suggest that the properties of GW190412 can be naturally explained by highergeneration mergers in an AGN disk. Due to the low mass ratio, we can expect that this may be a merger between a first-generation secondary BH with M 2 ≈ 8 M and a second-(or higher-) generation primary BH with M 1 ≈ 30 M . Indeed, q ∼ 0.28 and m bin ∼ 38 M is common for mergers in AGN disks (see Fig. 14 in Paper I). Furthermore, since mergers endow the remnant BH with high spin, the high value for the primary BH spin in GW190412 is consistent with it having experience one or more prior mergers. Finally, in our models, the AGN disk delivers BHs to the inner regions where binarysingle interactions frequently misalign the spins relative to the orbital angular momentum. If we include this event to the analysis in § 3.3.2, the Bayes factor between model M1 and the observed distribution for GW190412 is 0.70, which suggests that the AGN channel can naturally explain the properties of GW190412 well.

CONCLUSIONS
In this paper we investigated the distribution of the effective spin parameter χ eff for BH binaries merging in accretion disks of AGN. We performed one-dimensional N -body simulations, combined with semi-analytical prescriptions of the relevant processes. χ eff is enhanced by the alignment of BH spins toward the binary orbital angular momenta due to gas accretion, while it is reduced by the randomization of binary orbital angular momenta due to hard binary-single interactions. This is the first detailed estimate for the χ eff distribution of stellar-mass BH mergers in AGN disks, considering the effects of binary-single interactions and gas accretion. Our main results can be summarized as follows: 1. Due to the randomization of the binary orbital angular momentum directions by frequent binarysingle interactions, χ eff is symmetric around zero, if radial migration of BHs to the inner, densely populated regions is efficient. The median value of |χ eff | depends most strongly on the initial BH spin magnitudes and the efficiency of migration, and is much less impacted by the other parameters or prescriptions we considered.
2. The χ eff distribution for all observed events reported by the LIGO/Virgo collaborations and the IAS group during LIGO/Virgo O1 and O2 is roughly consistent with the distribution expected for mergers in AGN disks (Fig. 6). The KS probabilities between the χ eff distribution of all events and those in our models are typically ∼ 0.1 − 0.7 (P KS,χ eff , Table 2). The observed χ eff distribution slightly favors moderate values for the initial BH spins (|a 0 | 0.5; see Fig. 8).
3. Even for the worst-fitting model, the fractional contribution of mergers in AGN disks to all observed mergers is limited only to 0.3 (Fig. 9), and much higher contributions are allowed in our other models.
4. The positive correlation between m chirp and the dispersion of χ eff can be reproduced by AGNassisted mergers if the initial BH spin magnitude is low (Figs. 5 and 11, § 3.4). Also, mergers in AGN disks might be distinguished from other channels based on the chirp masses extending to values as high as ≈ 300 M (see also Paper I).
5. The properties of the recently announced gravitational-wave event, GW190412, including a low mass ratio, a high spin for the primary BH, and a spin component in the orbital plane, are naturally expected if it is a hierarchical merger in an AGN disk.
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No 638435 (GalNUC) and by the Hungarian National Research, Development, and Innovation Office grant NKFIH KH-125675. ZH acknowledges support from NASA grant NNX15AB19G and NSF grant 1715661. Simulations and analyses were carried out on Cray XC50 and computers at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. Table 2 The results in different models. The first two columns show the model number and indicate its variation from the fiducial model (M1). For example, in M16 ("No gas hard") binaries are not hardened by gas interaction, and M15 ("w/o gas mig.") excludes type I/II torques and the resulting radial migration in the AGN disk. In the next four columns, χ eff,med , χ eff,10 , χ eff,90 are the median, 10 percentile, and 90 percentile for the χ eff -and the median for the |χ eff |-distributions, respectively, in which observational errors are included. In the next three columns, P KS,χ eff , P KS,m chirp , and P KS,χ eff ,m chirp are, respectively, the KS probabilities that the all observed events were drawn from the χ eff -, m chirp -, and the joint (χ eff , m chirp )-distributions predicted in each model. In the next two columns, P KS,LV,χ eff and P KS,LV,m chirp are the KS probabilities that the events reported by the LIGO/Virgo collaborations were drawn from the predicted χ effand m chirp -distributions. In the last column, min i K A,obs,i is the lowest value of the Bayes factor among all observed GW events, evaluated for each event between the given model and the observed distribution itself.
input output Model Parameter χ eff,med χ eff,10 χ eff,90 |χ eff | med P KS,χ eff P KS,m chirp P KS,χ eff ,m chirp P KS,LV,χ eff P KS,LV,m chirp min i K   Zackay et al. (2019a). Note that reference [1] quotes the source-frame, whereas [2,3] quote the detector-frame chirp masses, together with their respective errors (columns 2 and 3, respectively). For the events found by [2-4], we calculate the dispersion of the source-frame chirp mass assuming no covariance between the parameters. K M1,obs,i is the Bayes factor between model M1 and the observed distribution for each event i.  log 10 (m chirp ) (t) observed dist. Figure 11. Same as Fig. 4, but the standard deviations are presented. Blue lines show the standard deviation of the mergers in each bin, and orange lines show the standard deviation of the χ eff distribution which is produced by performing a kernel-density estimate. The latter is shown as an estimate of the trend for the observed distribution in which the number of events is small.