# Mass-gap Mergers in Active Galactic Nuclei

Hiromichi Tagawa11affiliationmark:  , Bence Kocsis22affiliationmark:  , Zoltán
Haiman33affiliationmark:  , Imre Bartos44affiliationmark:  , Kazuyuki
Omukai11affiliationmark:  , Johan Samsing55affiliationmark:
11affiliationmark:  Astronomical Institute, Graduate School of Science, Tohoku
University, Aoba, Sendai 980-8578, Japan
22affiliationmark:  Rudolf Peierls Centre for Theoretical Physics, Clarendon
Laboratory, Parks Road, Oxford OX1 3PU, UK
33affiliationmark:  Department of Astronomy, Columbia University, 550 W. 120th
St., New York, NY, 10027, USA
44affiliationmark:  Department of Physics, University of Florida, PO Box
118440, Gainesville, FL 32611, USA
55affiliationmark:  Niels Bohr International Academy, The Niels Bohr
Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark  E-mail:
<EMAIL_ADDRESS>
###### Abstract

The recently discovered gravitational wave sources GW190521 and GW190814 have
shown evidence of BH mergers with masses and spins outside of the range
expected from isolated stellar evolution. These merging objects could have
undergone previous mergers. Such hierarchical mergers are predicted to be
frequent in active galactic nuclei (AGN) disks, where binaries form and evolve
efficiently by dynamical interactions and gaseous dissipation. Here we compare
the properties of these observed events to the theoretical models of mergers
in AGN disks, which are obtained by performing one-dimensional $N$-body
simulations combined with semi-analytical prescriptions. The high BH masses in
GW190521 are consistent with mergers of high-generation (high-g) BHs where the
initial progenitor stars had high metallicity, 2g BHs if the original
progenitors were metal-poor, or 1g BHs that had gained mass via super-
Eddington accretion. Other measured properties related to spin parameters in
GW190521 are also consistent with mergers in AGN disks. Furthermore, mergers
in the lower mass gap or those with low mass ratio as found in GW190814 and
GW190412 are also reproduced by mergers of 2g-1g or 1g-1g objects with
significant accretion in AGN disks. Finally, due to gas accretion, the massive
neutron star merger reported in GW190425 can be produced in an AGN disk.

###### Subject headings:

binaries: close – gravitational waves –galaxies: active – methods: numerical –
stars: black holes

##  1\. Introduction

Recently, several gravitational wave (GW) events were reported by LIGO and
Virgo whose measured physical properties pose interesting constraints on their
astrophysical origin. One of these is GW190521 (The LIGO Scientific
Collaboration et al., 2020b; Abbott et al., 2020d), in which the masses of the
merging BHs ($85^{+21}_{-14}\,{\rm M}_{\odot}$ and $66^{+17}_{-18}\,{\rm
M}_{\odot}$, but see Nitz & Capano 2020) indicate that the merging BHs fall
into the so-called “upper-mass gap”, as they are heavier than the maximum BH
mass imposed by (pulsation) pair-instability supernovae (Chatzopoulos &
Wheeler, 2012; Farmer et al., 2019). The components’ dimensionless spins
($a=0.69^{+0.27}_{-0.62}$ and $0.73^{+0.24}_{-0.64}$) are also larger than the
maximum natal spins predicted by stellar evolutionary models for angular
momentum transfer (Fuller & Ma, 2019). A possible electromagnetic counterpart
has been reported for GW190521, which, if real, would suggest mergers in gas-
rich environments (McKernan et al., 2019; Graham et al., 2020). Further
candidates for mergers in the upper-mass gap are reported in GW170729 (Abbott
et al., 2019), GW170817A (Zackay et al., 2019a), and six other candidates in
O3 (Abbott et al., 2020c). A possible explanation proposed for the large
masses in these events is isolated binary evolution, presuming that massive
BHs can form from Population III (Pop III) or metal-poor stars (Tanikawa et
al., 2020; Kinugawa et al., 2020; Farrell et al., 2020; Liu & Bromm, 2020; Liu
et al., 2020; Belczynski, 2020). However, due to their expected rapid
rotation, stellar remnants are unlikely to retain their hydrogen envelopes and
produce BHs as massive as in GW190521 (Chatzopoulos & Wheeler, 2012; Yoon et
al., 2012). Likewise, the probability that the progenitor of the $85\,{\rm
M}_{\odot}$ BH avoided a pulsation pair-instability supernova due to
uncertainties on the reaction rate of carbon burning is low (Farmer et al.,
2019). An alternative proposal is for Pop III stellar remnants to grow via
accretion (Safarzadeh & Haiman, 2020); this requires the remnant BHs to spend
a significant time in a dense gaseous environment of their host protogalaxies.

Similarly, the neutron star (NS) merger GW190425 indicates masses higher than
in NS binaries observed in the Milky Way (Abbott et al., 2020a). To explain
this event, a rapid merging pathway is proposed for the binary evolution
channel (Romero-Shaw et al., 2020; Safarzadeh et al., 2020b; Galaudage et al.,
2020). Additionally, mergers with very unequal masses have been reported –
GW190412 (The LIGO Scientific Collaboration et al., 2020a) and GW190814
(Abbott et al., 2020b) –, which are also atypical in stellar evolutionary
models of isolated binaries (Gerosa et al., 2020b; Olejak et al., 2020). The
mass of the secondary object of $2.6\,{\rm M}_{\odot}$ in GW190814 falls in
the so-called lower mass gap, predicted by stellar-evolutionary models (Fryer
et al. 2012; Zevin et al. 2020 but see Belczynski et al. 2012; Drozda et al.
2020) and obeyed by observations of X-ray binaries (Özel et al., 2012). The
object in the lower mass gap in GW190814 and a non-zero spin for the primary
BH ($0.43^{+0.16}_{-0.26}$) in GW190412 are consistent with a scenario in
which the merging compact objects (COs) had experienced previous mergers or
significant accretion. These events suggest that growth by gas accretion or
hierarchical mergers may be common among COs (see e.g. O’Leary et al., 2016;
Gerosa et al., 2020b; Abbott et al., 2020d; Hamers & Safarzadeh, 2020;
Rodriguez et al., 2020; Safarzadeh et al., 2020a; Safarzadeh & Haiman, 2020;
Yang et al., 2020a; Liu & Lai, 2020).

By analyzing the ensemble of events detected during LIGO/Virgo’s O1-O3a
observing runs, Kimball et al. (2020) and Tiwari & Fairhurst (2020) found
preference for at least one, but probably multiple hierarchical mergers in the
detected sample. The results by Kimball et al. (2020) strongly depend on the
escape velocity at the merger sites, with higher escape velocities favoring a
larger number of hierarchical mergers.

Active galactic nuclei (AGNs) represent environments where mass growth by gas
accretion is possible (Levin, 2007) and hierarchical mergers are frequent due
to efficient hardening by interaction with gas (e.g. Bartos et al., 2017;
Stone et al., 2017; McKernan et al., 2018, 2020b, 2020a), stars, and COs
(Tagawa et al., 2020b, hereafter Paper I). Additionally, the high escape
velocity in such systems helps retain recoiling merger remnants. Yang et al.
(2020b) investigated the possibility of lower mass gap mergers by assuming
that the objects are quickly delivered to migration traps of AGN disks where
they undergo mergers. In Paper I, we performed one-dimensional (1D) $N$-body
simulations combined with semi-analytical prescriptions for gaseous and
stellar effects, to follow the evolution of binary separation, along with the
radius and inclination of the orbit around the SMBH. The results showed that
binaries form efficiently during dynamical encounters in gas, and mergers are
facilitated by binary-single (BS) interactions in the gap-forming inner
regions. In Tagawa et al. (2020a, hereafter Paper II), we extended this model
to follow the evolution of BH spins, and binary orbital angular momentum
directions. In Tagawa et al. (2020c, hereafter Paper III), we also included
the effects of binary eccentricity, and found that mergers are often highly
eccentric in AGN disks if BS interactions are confined to the AGN plane due to
gaseous torques (see also Samsing et al. 2020) or if the objects migrate close
to the SMBH where they undergo GW capture events.

In this paper, we extend our previous models to include NSs. We then examine
the mass distribution of merging binaries in the AGN channel, and assess the
likelihood of mergers in the lower and upper mass gaps, those with highly
unequal mass ratios, and massive NS mergers as seen in GW190521, GW190814,
GW190412, and GW190425.

##  2\. Method

Our model is based on Papers I, II, and III, summarized in § 2.1, with new
ingredients described in $\S\,\ref{sec:new_objects}$.

###  2.1. Overview of model

We consider a model with a supermassive BH (SMBH) at the center of a galaxy,
surrounded by a gaseous accretion disk (hereafter AGN disk) and a spherical
nuclear star cluster (NSC), which are assumed to be in a fixed steady-state.
We follow the $N$-body evolution of COs consisting of those in the nuclear
stars cluster and those captured inside the AGN disk.

We follow several pathways to form binaries. Some fraction of COs are
initially in binaries. In the AGN disk, binaries form due to gas dynamical
friction during two-body encounters and dynamical interactions during three-
body encounters. Highly eccentric CO binaries also form due to the GW-capture
mechanism in single-single encounters (e.g. O’Leary et al., 2009) and BS
exchange interactions (Paper III).

Regarding the interaction with gas, the velocities of all COs relative to the
local AGN disk decrease due to accretion torque and gas dynamical friction.
The binaries’ semi-major axis decreases due to gas dynamical friction from the
AGN disk. Their semi-major axis and eccentricity evolve also due to type I/II
migration torque from a circumbinary disk. COs migrate toward the SMBH due to
type I/II torques from the AGN disk. CO masses, BH spins, and the orbital
angular momentum directions of binaries gradually evolve due gas accretion
(Paper II).

We also account for dynamical interactions with single stars/COs and CO
binaries. The binaries’ semi-major axis, velocities, orbital angular momentum
directions, and eccentricity evolve due to BS interactions as prescribed in
Leigh et al. (2018), and the velocities of all COs additionally evolve due to
scattering. The soft binaries, i.e. those with binding energy lower than that
of scattering objects, are either softened or ultimately disrupted during BS
interactions. Hard binaries become harder during BS interactions, and exchange
interactions may take place; the most massive pair is assumed to remain in the
binary after the interaction (Paper III).

The binaries’ separation and eccentricity decrease due to GW emission. When
the pericenter distance of the binary becomes smaller than the innermost
stable orbit, we assume that BHs merge, assign a kick velocity and mass loss
due to GW emission, and prescribe BH spin evolution at the merger as in Paper
II.

Table 1 The results of different models at $t=10\,\mathrm{Myr}$. The two leftmost (”input”) columns show the model number and the differences from the fiducial model. The ”output” columns list the fraction of mergers weighted by the detectable volume for BH-LGO ($f_{\mathrm{BHLGO}}$), LGO-LGO ($f_{\mathrm{LGOLGO}}$), LGO-NS ($f_{\mathrm{LGONS}}$), BH-BH ($f_{\mathrm{BHBH}}$), BH-NS ($f_{\mathrm{BHNS}}$), NS-NS ($f_{\mathrm{NSNS}}$) binaries, BH-BH binaries with $m_{\mathrm{bin}}>130\,{\rm M}_{\odot}$ ($f_{\mathrm{M130}}$), and the total number of mergers ($N_{\mathrm{mer}}$), where LGOs represent lower mass gap objects with masses in $2$–$5\,{\rm M}_{\odot}$.  input | output
---|---
Model | Parameter | $f_{\mathrm{BHLGO}}$ | $f_{\mathrm{LGOLGO}}$ | $f_{\mathrm{LGONS}}$ | $f_{\mathrm{BHBH}}$ | $f_{\mathrm{BHNS}}$ | $f_{\mathrm{NSNS}}$ | $f_{\mathrm{M130}}$ | $N_{\mathrm{mer}}$
M1 | Fiducial | $4\times 10^{-3}$ | $3\times 10^{-4}$ | $3\times 10^{-4}$ | 0.989 | $6\times 10^{-3}$ | $3\times 10^{-4}$ | $0.03$ | $1.6\times 10^{3}$
M2 | $m_{\mathrm{1gBH}}~{}\leq 45\,{\rm M}_{\odot}$ | $7\times 10^{-4}$ | $5\times 10^{-5}$ | $7\times 10^{-5}$ | 0.998 | $1\times 10^{-3}$ | $5\times 10^{-5}$ | 0.11 | $2.0\times 10^{3}$
M3 | $\Gamma_{\mathrm{Edd,cir}}=20$ | $5\times 10^{-3}$ | $5\times 10^{-5}$ | $9\times 10^{-5}$ | 0.994 | $8\times 10^{-4}$ | $5\times 10^{-5}$ | 0.35 | $1.9\times 10^{3}$

###  2.2. Evolution of neutron stars

While only stellar-mass BHs are followed in Papers I/II/III, here we
additionally follow the evolution of NSs among the $N$-body particles. We
assume that stars with zero-age main-sequence masses of $8$–$20\,{\rm
M}_{\odot}$ form NSs with $1.3\,{\rm M}_{\odot}$. There are initially $5\times
10^{4}$ NSs in the fiducial model. As NSs are lighter than BHs, we assume a
shallower radial density profile

$\displaystyle\frac{dN_{\rm NS,ini}(r)}{dr}\propto r^{-0.5},$ (1)

where $N_{\rm NS,ini}(r)$ is the total initial number of NSs within distance
$r$ from the central SMBH, and a higher velocity dispersion in the frame
comoving with the AGN ($\sigma_{\rm NS}=0.4v_{\mathrm{kep}}(r)$, where
$v_{\mathrm{kep}}(r)$ is the Keplerian velocity at $r$) compared to those for
BHs ($dN_{\rm BH,ini}(r)/dr\propto r^{0}$, $\sigma_{\rm
BH}=0.2v_{\mathrm{kep}}(r)$), as expected from relaxation processes (Hopman &
Alexander, 2006; Szolgyen & Kocsis, 2018). For simplicity, we assume that the
spin magnitudes of NSs are always zero, the mass lost during mergers is
$0.01\,{\rm M}_{\odot}$, and their radius is $10\,\mathrm{km}$, within which
NSs are assumed to merge in addition to the condition described above. We
regard objects with mass less than $2\,{\rm M}_{\odot}$ to be NSs, while we
call objects with masses in the range $2$–$5\,{\rm M}_{\odot}$ lower mass gap
objects (LGOs), and treat LGOs as BHs in the calculations. We assign no kicks
during NS-CO mergers for simplicity (but see Dietrich et al., 2017), while the
spin evolution during NS-CO mergers is prescribed as for BH-BH mergers but
with zero spin magnitudes for NSs. The interactions of NSs with gas, COs, and
stars are prescribed in the same manner as for BHs.

Figure 1.—  The parameter distributions for merging binaries weighted by the
detectable volume at $t=10\,\mathrm{Myr}$ in model M1. The four panels present
distributions of (a) the masses of primary and secondary object ($m_{1}$ vs.
$m_{2}$) in units of ${\rm M}_{\odot}$, (b) the dimensionless spins of primary
and secondary objects ($a_{1}$ vs. $a_{2}$), (c) the angles between the binary
orbital angular momentum and the BH spins ($\mathrm{cos}\theta_{1}$ and
$\mathrm{cos}\theta_{2}$), and (d) the mass-weighted sum of the individual
spin components perpendicular to the orbital plane ($\chi_{\mathrm{eff}}$) vs.
the precession parameter ($\chi_{\mathrm{p}}$), respectively. The error bars
correspond to the $90$ percentile credible intervals inferred in GW190521 (The
LIGO Scientific Collaboration et al., 2020b).  Figure 2.—  The parameter
distributions for merging binaries weighted by the detectable volume among
$n$g-$n$g mergers with $n\geq 4$ in model M1 (a), 2g-2g mergers in model M2
(b), and 1g-1g mergers in model M3 (c), where $n$g BHs represent BHs consisted
of merged $n$-1g BHs. Left, middle, and right panels show the distributions of
$m_{1}$ vs. $m_{2}$ in units of ${\rm M}_{\odot}$, and the dimensionless spin
parameters $a_{1}$ vs. $a_{2}$, and $\chi_{\mathrm{eff}}$ vs.
$\chi_{\mathrm{p}}$, respectively. The error bars correspond to the $90$
percentile credible intervals for GW190521 (The LIGO Scientific Collaboration
et al., 2020b).

###  2.3. Numerical choices

The parameters of the fiducial model (M1) are listed in Table 3 in the
Appendix, which are the same as those in Papers I/II/III. We also study the
additional models M2 and M3, which differ in the initial masses of stellar-
mass BHs and efficiency of gas accretion.

In model M1, we set the mass of first-generation (1g) BHs to
$m_{\mathrm{1gBH}}\leq 15\,{\rm M}_{\odot}$ given the observed high
metallicity population for $90\%$ of stars in the Milky Way NSC (e.g. Do et
al., 2020) and the observations indicating strong mass loss during the
luminous blue variable phase of solar-metallicity stars (e.g. Chen et al.,
2015), which likely limits the maximum mass of 1g BHs to be $\sim 15\,{\rm
M}_{\odot}$ (e.g. Belczynski et al., 2010).

On the other hand, in model M2, the initial masses of 1g BHs are multiplied by
3 compared to the values in model M1, and so $m_{\mathrm{1gBH}}\leq 45\,{\rm
M}_{\odot}$. Such high masses for 1g BHs may exist as indicated by a
moderately metal-poor subpopulation in the Milky Way’s NSC (Do et al., 2020),
and given that there are significant uncertainties in the initial mass
function of 1g BHs in NSCs due to star formation in strong gravitational field
(e.g. Nayakshin et al., 2007), gas accretion (e.g. Davies & Lin, 2020),
uncertainties in stellar evolution (e.g. Belczynski et al., 2010), and the
diversity in the possible origin of objects (e.g. Mapelli & Gualandris, 2016).

In model M3, we increase the maximum accretion rate in Eddington units onto
stellar-mass COs to $20$, given the uncertainties due to radiation feedback
onto dust (Toyouchi et al., 2020) and rotating gas (Jiang et al., 2014;
Sugimura et al., 2018), and the possible anisotropy of radiation (Proga et
al., 2000).

Figure 3.—  The distribution of $m_{1}$ and $m_{2}$ in units of ${\rm
M}_{\odot}$ weighted by the detectable volume in models M1–M3 at
$t=100\,\mathrm{Myr}$. The error bars correspond to the $90$ percentile
credible intervals for GW190814 (panels a–c, Abbott et al. 2020b), GW190412
(panels d and e, The LIGO Scientific Collaboration et al. 2020a), and GW190425
(panel f, Abbott et al. 2020a) for low-(red) and high-(orange) spin priors.

##  3\. Results

###  3.1. Upper mass-gap mergers

We first compare the properties of predicted mergers in AGN disks with those
of the upper-mass gap event GW190521. The parameter distributions for all
events at $t=10\,\mathrm{Myr}$ in model M1 are shown in Fig. 1, while those
for specific generation mergers at $t=100\,\mathrm{Myr}$ reproducing well
GW190521 in models M1–M3 in Fig. 2. Note that the maximum binary masses are
not sensitive to time in later phases of $t\geq 10\,\mathrm{Myr}$ as discussed
in the Appendix. In Fig. 2, results at $t=100\,\mathrm{Myr}$ are presented to
reduce variation due to a small number of mergers. The distributions are
weighted by the observable volume depending on the masses of merging binaries
calibrated for aLIGO Handford prior to O3 as in Paper I, while the dependence
on either $\chi_{\mathrm{eff}}$ or eccentricity is ignored in terms of
detection volume for simplicity. The plots are smoothed by performing kernel-
density estimate using Gaussian kernels whose bandwidth is chosen to satisfy
Scott’s rule (Scott, 1992).

In model M1, the primary-mass distribution extends to as high as $>100{\rm
M}_{\odot}$ despite the low initial BH mass of $\leq 15\,{\rm M}_{\odot}$
(Fig. 1 a). High-mass events with $\gtrsim 50~{}{\rm M_{\odot}}$ BHs,
including GW190521, primarily arise from high-g mergers as seen in Fig. 2 a,
where the distribution is shown for the mergers between high $n$-generation
objects with $n\geq 4$. In model M1 at $t=10\,\mathrm{Myr}$, the expected
detection fraction of upper-mass gap mergers with $50\,{\rm M}_{\odot}\leq
m_{1}\leq 150\,{\rm M}_{\odot}$ contributed by $3$, 4, 5, 6, 7, 8, and $n\geq
9$th-generation primary BHs is 0, 5, 17, 30, 21, 3, and 24 percent,
respectively, where an $n$g BH signifies a BH consisted of $n$ 1g BHs as the
merged progenitors. As the mass of $n$g BHs is at most $\sim n\times 15\,{\rm
M}_{\odot}$ in model M1, primary BHs of GW190521-like mergers $m_{1}\sim
85\,{\rm M}_{\odot}$ come from $n\geq 6$-g BHs. The fraction of massive binary
mergers with $m_{\mathrm{bin}}\geq 130\,{\rm M}_{\odot}$ ($\sim 90$ percentile
lower limit for GW190521) over all events ($f_{\mathrm{M130}}$) is 0.03 in
model M1. Since GW190521 was found among $\sim 60$ triggers, the rate of this
type of events can be explained by the AGN channel if it is the dominant
channel for all the detected events. Here, note that the maximum mass of
merging binaries ($\sim 185\,{\rm M}_{\odot}$) is significantly reduced from
Papers I ($\sim 10^{4}\,{\rm M}_{\odot}$), which is discussed in the Appendix.

The mergers tend to have high spin around $\sim 0.7$ (Fig. 1 b). This is
because the spin of their members is already high as a result of previous
episode of mergers (e.g. Buonanno et al., 2008). Also, the angles
$\theta_{1}$, $\theta_{2}$ between the binary orbital angular momentum and the
primary and secondary BH spin, respectively, are widely distributed in the
entire range of ($0$–$\pi$) (Fig. 1 c) due to several BS interactions, which
are assumed to randomize the orientation of the orbital angular momentum. Yet,
strong correlation between $\theta_{1}$ and $\theta_{2}$ is still observed
owing to spin alignment by gas accretion. High-g mergers have a wide range of
$\chi_{\mathrm{p}}$ and $\chi_{\mathrm{eff}}$ (right panel of Fig. 2 a) due to
the spin enhancement during mergers and the randomization of the orbital
angular momentum directions by BS interactions.

In summary, the results above demonstrate that mergers between high-g BHs in
model M1 is consistent with the observed properties ($\chi_{\mathrm{eff}}$,
$\chi_{\mathrm{p}}$, $m_{1}$, $m_{2}$, $a_{1}$, $a_{2}$, $\theta_{1}$,
$\theta_{2}$, rate) of GW190521.

Fig. 2b shows the parameter distributions of the 2g-2g BH mergers in model M2,
in which the maximum initial BH mass is elevated to $45\,{\rm M}_{\odot}$,
typical for metal-poor progenitors. This high initial mass allows 2g-2g BH
mergers to reach the mass range of GW190521. The fraction of mergers with
$m_{\mathrm{bin}}\geq 130\,{\rm M}_{\odot}$ becomes as high as 0.11 in this
model. Since some stars in NSCs are indeed metal-poor (Do et al., 2020), the
observed rate of GW190521-like events is in concordance also with model M2 as
long as the AGN channel comprises at least $\sim 15\%$ of all the GW events
detected to date. This number is calculated from the observed massive merger
fraction ($\sim 1/60$) divided by the fraction in our model ($0.11$).

Fig. 2c shows the results for the 1g-1g mergers in model M3, in which the
accretion rate can be as high as 20 times the Eddington accretion rate. Now,
owing to rapid growth by accretion, even some 1g-1g mergers reach the mass
range of GW190521. The gas accretion also significantly enhances the BH spin,
and $a_{1}$, $a_{2}$, and $\chi_{\mathrm{p}}$ (or $|\chi_{\mathrm{eff}}|$).
Such rapid accretion scenario may be favored by the presence of events with
high $|\chi_{\mathrm{eff}}|$ values such as
$\chi_{\mathrm{eff}}=0.81^{+0.15}_{-0.21}$ of GW151216 (Zackay et al., 2019b).
Precise estimate for the BH spins in future will further constrain its
validity.

In conclusion, the binary property, including the mass, spin and occurence
rate, of GW190521 is consistent with the events via the AGN channel. In
particular, its origin can be well explained either by high($\geq 4$)-g merger
of massive BHs, by 2g BH mergers which have massive low-metallicity progenitor
stars, or by mergers of BHs that have rapidly grown by super-Eddington
accretion.

###  3.2. lower mass-gap, low mass-ratio, and massive NS mergers

Next, we examine whether our model can reproduce lower mass-gap, low mass-
ratio, and massive NS merger events, like GW190814, GW190412, and GW190425,
respectively. Fig. 3 shows the mass distribution of the merging objects for
models M1–M3. The mass ranges shown in each panel is chosen so that the
distribution around the corresponding observed events is clearly seen. We do
not show the result of model M2 for GW190412 since the initial BH mass falls
in the range $15$–$45\,{\rm M}_{\odot}$ by construction and thus is in
contradiction with the secondary mass of the event
($m_{2}=8.4^{+1.7}_{-1.0}\,{\rm M}_{\odot}$).

In models M1–M3, moderate probability ($f_{\mathrm{BHLGO}}=0.07$–$0.4\%$,
Table 1) of the detection fraction is predicted for mergers between LGOs and
BHs. Note that $f_{\mathrm{BHLGO}}$ can be as high as $\gtrsim 1\%$ depending
on the adopted parameters (see Appendix). This is marginally consistent with
the rate of events like GW190814, which was found from $\sim 60$ triggers in
LIGO/Virgo observations, provided the AGN channel is the dominant channel for
all the detected GW events. Secondary objects of around $2.6\,{\rm M}_{\odot}$
can be produced by mergers of NSs (panels a and b), while in model M3, 1g
objects can be as massive as this thanks to super-Eddington accretion. With
low spin magnitude ($\lesssim 0.07$, Abbott et al. 2020b), the primary BH of
GW190814 is probably 1g and has hardly grown by accretion. This fact favors a
merger between a massive-1g BH and a LGO as the origin of GW190814 as in model
M2.

Mergers like GW190412, having the low mass-ratio ($q=0.28^{+0.13}_{-0.07}$)
and the high primary-BH spin ($a_{1}=0.43^{+0.16}_{-0.26}$), can be explained
by $2$g-$1$g mergers in model M1 (panel d of Fig. 3) and by 1g-1g mergers in
model M3 (panel e). The detection probability of $2$g-$1$g mergers is $\sim
16\%$ in model M1, while that for low mass-ratio mergers with $q<0.3$ is $\sim
8\%$ in model M3. Also, non-zero spin for the primary BH in GW190412 may be
consistent with being $2$g (or high-g) as in model M1, or experienced
significant accretion as in model M3.

Overall, the properties of GW190814/GW190412 can be explained by the AGN
channel either by mergers among $2$g-$1$g objects or mergers of BHs grown by
super-Eddington accretion. The possible 1g-$20\,{\rm M}_{\odot}$ BH of
GW190814 prefers moderately metal-poor progenitors. GW190412 allows for solar-
metallicity progenitors because high mass of the primary BH ($\sim 30\,{\rm
M}_{\odot}$) suggests it to be a 2 or higher-g object or growth by significant
accretion. As both high- and low-metallicity stars are present in NSCs (Do et
al., 2020), AGN disks can provide candidate sites for both events.

Finally, the massive NS merger event GW190425 is well reproduced by growth of
the progenitors by gas accretion in model M1 at $t=100\,\mathrm{Myr}$ (panel f
of Fig. 3). This result is consistent with the finding in Yang et al. (2020a)
that gas accretion expected in the AGN channel can produce GW190425-like
massive NS mergers. In model M1, the merging NS mass extends to $\sim
1.6\,{\rm M}_{\odot}$ at $t=10\,\mathrm{Myr}$, slightly lower than the CO
masses of GW190425. Since the mass accreted in a fixed time interval is
roughly proportional to the gas accretion rate, which is here set to the
Eddington rate with $10\%$ radiative efficiency, the mass distribution of NS
mergers can constrain the typical accretion rate onto the progenitor and the
AGN disk lifetime once a large number of NS mergers in AGNs are observed.
However, the detection fraction of NS-NS mergers over all mergers is
$f_{\mathrm{NSNS}}\sim 5\times 10^{-5}$–$3\times 10^{-3}$ (Table 1, see also
Appendix), which is much lower than the actual observed fraction ($\sim
1/30$). Therefore, it might be unlikely that GW190425 is originated from AGNs.

We have also calculated several models with different sets of the parameters
related to the properties of the AGN disk, BH population, etc. and found that
the fraction of mass-gap mergers is only moderately affected by the difference
in their choice (see Appendix).

###  3.3. Influence of assumptions

Here, we discuss the main assumptions in our model and their possible
influence on our main conclusions. First, we adopted the steady-state AGN disk
model proposed by Thompson et al. (2005) and the initial three-dimensional
number density and mass distribution for COs as inferred for the Galactic
Center (e.g. Feldmeier et al., 2014; Schödel et al., 2014; Do et al., 2020).
Since there are uncertainties in the properties of the AGN disks and the CO
distributions, and they are also expected to vary from one galaxy to another,
we investigated the dependence of the results on the gas accretion rate from
outer radius of the simulation (models M11–M13), the size of the AGN disks
(models M4 and M10), the mass of the central SMBH (model M14), the initial
mass function of COs (model M15), the stellar mass enclosed within
$3\,\mathrm{pc}$ (model M16), the initial density profile of COs (model M17),
and the initial velocity dispersion of COs (models M5, M18, M19) in the
Appendix.

Although the occurrence rates of specific events like mergers moderately
depend on these parameters (Table 3 in the Appendix), the range of variations
does not significantly affect our conclusions that frequent hierarchical
mergers in the AGN disks are consistent with properties of the observed merger
events. On the other hand, we have not modelled the time evolution of the
structure of the AGN disk, the formation and evolution of compact objects (and
their binaries) other than BHs and NSs, and the possible presence of massive
perturbers (Deme et al., 2020), which may affect various properties of
mergers.

Second, we assumed that COs embedded in the AGN disk migrate toward the SMBH
following the formulae for Type I/II migration. However, it is not clear
whether radial migration operates as prescribed by these planetary formulae,
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). To investigate the impact of the migration on the
results, we performed model M8, in which the radial migration is assumed to be
inefficient. In this model, the number and the maximum mass of mergers are
reduced as the rate of binary-single interactions, which facilitate mergers,
is strongly enhanced by the migration to the inner regions of the AGN disk,
where gaps are opened, and the space density of BHs is high (see Paper I).

Third, hierarchical triples composed of three stellar-mass compact objects
(e.g. Antonini et al. 2017; see also Section 5.7.2 of Paper I), neglected in
this study, may influence the spin distribution for some fraction of mergers
and the merging mass distribution.

Fourth, we employed a semi-analytical prescription for binary-single
interactions. However, this prescription may not be a good approximation for
low mass-ratio interactions. Also, we assumed that the most massive pair
remains in a binary after hard binary-single interactions, and softer binaries
are always disrupted after hard binary-binary interactions for simplicity. By
calculating the interactions using $N$-body simulations, the merging mass
distribution may be significantly modified.

As impacts of several of these assumptions on results are uncertain, these
should be taken into account in future works.

##  4\. Conclusions

In this paper, we have followed the evolution of compact objects in AGN disks
by way of one-dimensional $N$-body simulations, combined with semi-analytical
prescriptions incorporating gaseous and stellar interactions, by extending our
previous model to include not only black holes but also neutron stars. We have
then compared the expected properties of compact-object mergers in AGN disks
with those of the observed events GW190521, GW190814, GW190412, and GW190425.

Our main findings are summarized as follows:

  1. 1.

The AGN channel can produce high BH-mass mergers like GW190521 in multiple
ways: (i) mergers among high-g BHs, (ii) mergers of 2g-2g BHs with metal-poor
progenitors, or (iii) mergers of BHs which have acquired mass via moderately
super-Eddington accretion. Events with roughly equal-mass ratio are expected
in this channel, unlike those occurring in the migration traps of the AGN
disks.

  2. 2.

Mergers in the lower mass-gap or with unequal masses like GW190814 and
GW190412 can be explained in the AGN channel by 2g-1g mergers or 1g-1g mergers
with significant gas accretion, with moderate occurrence probability.

  3. 3.

The Massive neutron-star merger GW190425 can be explained by the AGN channel
if the progenitors had acquired significant mass via accretion, although the
predicted rate is much lower than the observed rate.

The authors thank Will Farr for useful suggestions. This work is financially
supported by the Grants-in-Aid for Basic Research by the Ministry of
Education, Science and Culture of Japan (HT:17H01102, 17H06360, KO:17H02869,
17H01102, 17H06360). ZH acknowledges support from NASA grant NNX15AB19G and
NSF grants 1715661 and 2006176. This work received founding from the European
Research Council (ERC) under the European Union’s Horizon 2020 Programme for
Research and Innovation ERC-2014-STG under grant agreement No. 638435 (GalNUC)
(to BK). J.S. is supported by the European Unions Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie grant agreement No.
844629. IB acknowledges support from the Alfred P. Sloan Foundation and from
the University of Florida. Simulations and analyses were carried out on Cray
XC50 and computers at the Center for Computational Astrophysics, National
Astronomical Observatory of Japan.

## References

  * Abbott et al. (2019) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Physical Review X, 9, 031040
  * Abbott et al. (2020a) —. 2020a, ApJ, 892, L3
  * Abbott et al. (2020b) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b, ApJ, 896, L44
  * Abbott et al. (2020c) —. 2020c, arXiv e-prints, arXiv:2010.14527
  * Abbott et al. (2020d) —. 2020d, ApJ, 900, L13
  * Antonini et al. (2017) Antonini, F., Toonen, S., & Hamers, A. S. 2017, ApJ, 841, 77
  * Bartos et al. (2017) Bartos, I., Kocsis, B., Haiman, Z., & Márka, S. 2017, ApJ, 835, 165
  * Baruteau & Lin (2010) Baruteau, C., & Lin, D. N. C. 2010, ApJ, 709, 759
  * Belczynski (2020) Belczynski, K. 2020, arXiv e-prints, arXiv:2009.13526
  * Belczynski et al. (2010) Belczynski, K., Bulik, T., Fryer, C., et al. 2010, ApJ, 714, 1217
  * Belczynski et al. (2012) Belczynski, K., Wiktorowicz, G., Fryer, C. L., Holz, D. E., & Kalogera, V. 2012, ApJ, 757, 91
  * Belczynski et al. (2018) Belczynski, K., Bulik, T., Olejak, A., et al. 2018, arXiv e-prints, arXiv:1812.10065
  * Bellovary et al. (2016) Bellovary, J. M., Mac Low, M. M., McKernan, B., & Ford, K. E. S. 2016, ApJ, 819, L17
  * Broz et al. (2018) Broz, M., Chrenko, O., Nesvornu, D., & Lambrechts, M. 2018, A&A, 620, A157
  * Buonanno et al. (2008) Buonanno, A., Kidder, L. E., & Lehner, L. 2008, Phys. Rev. D, 77, 026004
  * Chatzopoulos & Wheeler (2012) Chatzopoulos, E., & Wheeler, J. C. 2012, ApJ, 748, 42
  * Chen et al. (2015) Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068
  * Czerny et al. (2019) Czerny, B., Wang, J.-M., Du, P., et al. 2019, ApJ, 870, 84
  * Davies & Lin (2020) Davies, M. B., & Lin, D. N. C. 2020, MNRAS, arXiv:2008.10033 [astro-ph.GA]
  * del Valle & Volonteri (2018) del Valle, L., & Volonteri, M. 2018, MNRAS, 480, 439
  * Deme et al. (2020) Deme, B., Meiron, Y., & Kocsis, B. 2020, ApJ, 892, 130
  * Dietrich et al. (2017) Dietrich, T., Ujevic, M., Tichy, W., Bernuzzi, S., & Brügmann, B. 2017, Phys. Rev. D, 95, 024029
  * Do et al. (2020) Do, T., David Martinez, G., Kerzendorf, W., et al. 2020, ApJ, 901, L28
  * Drozda et al. (2020) Drozda, P., Belczynski, K., O’Shaughnessy, R., Bulik, T., & Fryer, C. L. 2020, arXiv e-prints, arXiv:2009.06655
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53
  * Farrell et al. (2020) Farrell, E. J., Groh, J. H., Hirschi, R., et al. 2020, arXiv e-prints, arXiv:2009.06585
  * Feldmeier et al. (2014) Feldmeier, A., Neumayer, N., Seth, A., et al. 2014, A&A, 570, A2
  * Fragione et al. (2019) Fragione, G., Grishin, E., Leigh, N. W. C., Perets, H. B., & Perna, R. 2019, MNRAS, 488, 47
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91
  * Fuller & Ma (2019) Fuller, J., & Ma, L. 2019, ApJ, 881, L1
  * Galaudage et al. (2020) Galaudage, S., Adamcewicz, C., Zhu, X.-J., & Thrane, E. 2020, arXiv e-prints, arXiv:2011.01495
  * Gayathri et al. (2020) Gayathri, V., Healy, J., Lange, J., et al. 2020, arXiv e-prints, arXiv:2009.05461
  * Gerosa et al. (2020a) Gerosa, D., Mould, M., Gangardt, D., et al. 2020a, arXiv e-prints, arXiv:2011.11948
  * Gerosa et al. (2020b) Gerosa, D., Vitale, S., & Berti, E. 2020b, arXiv e-prints, arXiv:2005.04243
  * Gondán & Kocsis (2020) Gondán, L., & Kocsis, B. 2020, arXiv e-prints, arXiv:2011.02507
  * Graham et al. (2020) Graham, M. J., Ford, K. E. S., McKernan, B., et al. 2020, Phys. Rev. Lett., 124, 251102
  * Greene & Ho (2007) Greene, J. E., & Ho, L. C. 2007, ApJ, 667, 131
  * Hamers & Safarzadeh (2020) Hamers, A. S., & Safarzadeh, M. 2020, ApJ, 898, 99
  * Hopkins et al. (2005) Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2005, ApJ, 630, 705
  * Hopman & Alexander (2006) Hopman, C., & Alexander, T. 2006, ApJL, 645, L133
  * Jiang et al. (2014) Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2014, ApJ, 796, 106
  * Kimball et al. (2020) Kimball, C., Talbot, C., Berry, C. P. L., et al. 2020, arXiv e-prints, arXiv:2011.05332
  * Kinugawa et al. (2020) Kinugawa, T., Nakamura, T., & Nakano, H. 2020, arXiv e-prints, arXiv:2009.06922
  * Laughlin et al. (2004) Laughlin, G., Steinacker, A., & Adams, F. C. 2004, ApJ, 608, 489
  * Leigh et al. (2018) Leigh, N. W. C., Geller, A. M., McKernan, B., et al. 2018, MNRAS, 474, 5672
  * Levin (2007) Levin, Y. 2007, MNRAS, 374, 515
  * Liu & Bromm (2020) Liu, B., & Bromm, V. 2020, ApJ, 903, L40
  * Liu & Lai (2020) Liu, B., & Lai, D. 2020, arXiv e-prints, arXiv:2009.10068
  * Liu et al. (2020) Liu, B., Meynet, G., & Bromm, V. 2020, MNRAS, arXiv:2010.05824 [astro-ph.GA]
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415
  * Mapelli & Gualandris (2016) Mapelli, M., & Gualandris, A. 2016, Star Formation and Dynamics in the Galactic Centre, Vol. 905 (Lecture Notes in Physics, Springer, Cham), 205
  * McKernan et al. (2020a) McKernan, B., Ford, K. E. S., & O’Shaughnessy, R. 2020a, MNRAS, 498, 4088
  * McKernan et al. (2020b) McKernan, B., Ford, K. E. S., O’Shaugnessy, R., & Wysocki, D. 2020b, MNRAS, 494, 1203
  * McKernan et al. (2018) McKernan, B., Ford, K. E. S., Bellovary, J., et al. 2018, ApJ, 866, 66
  * McKernan et al. (2019) McKernan, B., Ford, K. E. S., Bartos, I., et al. 2019, ApJ, 884, L50
  * Nayakshin et al. (2007) Nayakshin, S., Cuadra, J., & Springel, V. 2007, MNRAS, 379, 21
  * Nitz & Capano (2020) Nitz, A. H., & Capano, C. D. 2020, arXiv e-prints, arXiv:2010.12558
  * O’Leary et al. (2009) O’Leary, R. M., Kocsis, B., & Loeb, A. 2009, MNRAS, 395, 2127
  * O’Leary et al. (2016) O’Leary, R. M., Meiron, Y., & Kocsis, B. 2016, ApJ, 824, L12
  * Olejak et al. (2020) Olejak, A., Fishbach, M., Belczynski, K., et al. 2020, ApJ, 901, L39
  * Özel et al. (2012) Özel, F., Psaltis, D., Narayan, R., & Santos Villarreal, A. 2012, ApJ, 757, 55
  * Proga et al. (2000) Proga, D., Stone, J. M., & Kallman, T. R. 2000, ApJ, 543, 686
  * Regan et al. (2019) Regan, J. A., Downes, T. P., Volonteri, M., et al. 2019, arXiv e-prints, arXiv:1811.04953
  * Rodriguez et al. (2020) Rodriguez, C. L., Kremer, K., Grudić, M. Y., et al. 2020, ApJ, 896, L10
  * Romero-Shaw et al. (2020) Romero-Shaw, I. M., Farrow, N., Stevenson, S., Thrane, E., & Zhu, X.-J. 2020, MNRAS, 496, L64
  * Safarzadeh & Haiman (2020) Safarzadeh, M., & Haiman, Z. 2020, arXiv e-prints, arXiv:2009.09320
  * Safarzadeh et al. (2020a) Safarzadeh, M., Hamers, A. S., Loeb, A., & Berger, E. 2020a, ApJ, 888, L3
  * Safarzadeh et al. (2020b) Safarzadeh, M., Ramirez-Ruiz, E., & Berger, E. 2020b, ApJ, 900, 13
  * Samsing et al. (2020) Samsing, J., Bartos, I., D’Orazio, D. J., et al. 2020, arXiv e-prints, arXiv:2010.09765
  * Schödel et al. (2014) Schödel, R., Feldmeier, A., Kunneriath, D., et al. 2014, A&A, 566, A47
  * Scott (1992) Scott, D. 1992, Multivariate Density Estimation: Theory, Practice, and Visualization, A Wiley-interscience publication (Wiley)
  * Shankar et al. (2013) Shankar, F., Weinberg, D. H., & Miralda-Escudé, J. 2013, MNRAS, 428, 421
  * Stone et al. (2017) Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, MNRAS, 464, 946
  * Sugimura et al. (2018) Sugimura, K., Hosokawa, T., Yajima, H., Inayoshi, K., & Omukai, K. 2018, MNRAS, 478, 3961
  * Szolgyen & Kocsis (2018) Szolgyen, A., & Kocsis, B. 2018, Phys. Rev. Lett., 121, 101101
  * Tagawa et al. (2020a) Tagawa, H., Haiman, Z., Bartos, I., & Kocsis, B. 2020a, ApJ, 899, 26
  * Tagawa et al. (2020b) Tagawa, H., Haiman, Z., & Kocsis, B. 2020b, ApJ, 898, 25
  * Tagawa et al. (2020c) Tagawa, H., Kocsis, B., Haiman, Z., et al. 2020c, arXiv e-prints, arXiv:2010.10526
  * Tanikawa et al. (2020) Tanikawa, A., Kinugawa, T., Yoshida, T., Hijikawa, K., & Umeda, H. 2020, arXiv e-prints, arXiv:2010.07616
  * The LIGO Scientific Collaboration et al. (2020a) The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, R., et al. 2020a, arXiv e-prints, arXiv:2004.08342
  * The LIGO Scientific Collaboration et al. (2020b) —. 2020b, arXiv e-prints, arXiv:2009.01075
  * The LIGO Scientific Collaboration et al. (2020c) —. 2020c, arXiv e-prints, arXiv:2010.14533
  * Thompson et al. (2005) Thompson, T. A., Quataert, E., & Murray, N. 2005, ApJ, 630, 167
  * Tiwari & Fairhurst (2020) Tiwari, V., & Fairhurst, S. 2020, arXiv e-prints, arXiv:2011.04502
  * Toyouchi et al. (2020) Toyouchi, D., Hosokawa, T., Sugimura, K., & Kuiper, R. 2020, MNRAS, 496, 1909
  * Yang et al. (2019) Yang, Y., Bartos, I., Haiman, Z., et al. 2019, ApJ, 876, 122
  * Yang et al. (2020a) Yang, Y., Gayathri, V., Bartos, I., et al. 2020a, arXiv e-prints, arXiv:2007.04781
  * Yang et al. (2020b) —. 2020b, arXiv e-prints, arXiv:2007.04781
  * Yoon et al. (2012) Yoon, S. C., Dierks, A., & Langer, N. 2012, A&A, 542, A113
  * Zackay et al. (2019a) Zackay, B., Dai, L., Venumadhav, T., Roulet, J., & Zaldarriaga, M. 2019a, arXiv e-prints, arXiv:1910.09528
  * Zackay et al. (2019b) Zackay, B., Venumadhav, T., Dai, L., Roulet, J., & Zaldarriaga, M. 2019b, Phys. Rev. D, 100, 023007
  * Zevin et al. (2020) Zevin, M., Spera, M., Berry, C. P. L., & Kalogera, V. 2020, ApJ, 899, L1

Table 2 Same as Table 1, but results for all models. “No gas hardening” and “No gas migration” respectively represent the models in which the gaseous torques are neglected with respect to the evolution of the binary semi-major axis and the radial position of the binary center of mass within the AGN disk. “BS in 2D” and “BS in 2$\rightarrow$3D” represents the models in which BS interactions occur in 2D and 2$\rightarrow$3D space, respectively (Paper III).  input | output
---|---
Model | Parameter | $f_{\mathrm{BHLGO}}$ | $f_{\mathrm{LGOLGO}}$ | $f_{\mathrm{LGONS}}$ | $f_{\mathrm{BHBH}}$ | $f_{\mathrm{BHNS}}$ | $f_{\mathrm{NSNS}}$ | $f_{\mathrm{M130}}$ | $N_{\mathrm{mer}}$
M1 | Fiducial | $4\times 10^{-3}$ | $3\times 10^{-4}$ | $3\times 10^{-4}$ | 0.989 | $6\times 10^{-3}$ | $3\times 10^{-4}$ | $0.03$ | $1.6\times 10^{3}$
M2 | $m_{\mathrm{1gBH}}~{}\times=3$ | $7\times 10^{-4}$ | $5\times 10^{-5}$ | $7\times 10^{-5}$ | 0.998 | $1\times 10^{-3}$ | $5\times 10^{-5}$ | 0.11 | $2.0\times 10^{3}$
M3 | $\Gamma_{\mathrm{Edd,cir}}=20$ | $5\times 10^{-3}$ | $5\times 10^{-5}$ | $9\times 10^{-5}$ | 0.994 | $8\times 10^{-4}$ | $5\times 10^{-5}$ | 0.35 | $1.9\times 10^{3}$
M4 | $r_{\mathrm{out,BH}}=0.03\,\mathrm{pc}$ | $0.10$ | $2\times 10^{-3}$ | $3\times 10^{-3}$ | 0.82 | $0.07$ | $3\times 10^{-3}$ | 0 | $3.1\times 10^{2}$
M5 | ${\beta}_{\mathrm{v,BH}}={\beta}_{\mathrm{v,NS}}=1$ | $0.03$ | $2\times 10^{-3}$ | $2\times 10^{-3}$ | 0.94 | $0.03$ | $2\times 10^{-3}$ | 0 | $5.9\times 10^{2}$
M6 | $f_{\mathrm{LGO}}=0.1$ | $0.01$ | $1\times 10^{-3}$ | $7\times 10^{-4}$ | 0.984 | $4\times 10^{-3}$ | $4\times 10^{-4}$ | 0.03 | $1.8\times 10^{3}$
M7 | $f_{\mathrm{LGO}}=0.01$ | $6\times 10^{-3}$ | $3\times 10^{-4}$ | $5\times 10^{-4}$ | 0.989 | $5\times 10^{-3}$ | $4\times 10^{-4}$ | 0.01 | $1.6\times 10^{3}$
M8 | No gas mig. | $0.02$ | $1\times 10^{-3}$ | $2\times 10^{-3}$ | 0.95 | $0.03$ | $2\times 10^{-3}$ | 0 | $7.1\times 10^{2}$
M9 | No gas hard. | $2\times 10^{-3}$ | $7\times 10^{-4}$ | $2\times 10^{-3}$ | 0.98 | $0.01$ | $1\times 10^{-3}$ | 0.11 | $1.5\times 10^{3}$
M10 | $r_{\mathrm{out,BH}}=0.3\,\mathrm{pc}$ | $8\times 10^{-3}$ | $7\times 10^{-4}$ | $1\times 10^{-3}$ | 0.98 | $0.01$ | $8\times 10^{-4}$ | 0.09 | $1.0\times 10^{3}$
M11 | ${\dot{M}}_{\mathrm{out}}={\dot{M}}_{\mathrm{Edd}}$ | $6\times 10^{-3}$ | $2\times 10^{-4}$ | $5\times 10^{-4}$ | 0.988 | $5\times 10^{-3}$ | $4\times 10^{-4}$ | 0.06 | $1.7\times 10^{3}$
M12 | ${\dot{M}}_{\mathrm{out}}=0.01{\dot{M}}_{\mathrm{Edd}}$ | $2\times 10^{-3}$ | $5\times 10^{-4}$ | $6\times 10^{-4}$ | 0.992 | $4\times 10^{-3}$ | $7\times 10^{-4}$ | 0.01 | $8.9\times 10^{2}$
M13 | ${\dot{M}}_{\mathrm{out}}=0.001{\dot{M}}_{\mathrm{Edd}}$ | $3\times 10^{-3}$ | $0$ | $4\times 10^{-4}$ | 0.98 | $0.02$ | $1\times 10^{-3}$ | 0 | $3.1\times 10^{2}$
M14 | $M_{\mathrm{SMBH}}=4\times 10^{7}{\rm M}_{\odot}$ | $8\times 10^{-3}$ | $2\times 10^{-4}$ | $4\times 10^{-4}$ | 0.98 | $7\times 10^{-3}$ | $3\times 10^{-4}$ | 0.07 | $1.6\times 10^{3}$
M15 | $\delta_{\mathrm{IMF}}=1.7$ | $1\times 10^{-3}$ | $1\times 10^{-4}$ | $1\times 10^{-4}$ | 0.997 | $1\times 10^{-3}$ | $1\times 10^{-4}$ | 0.09 | $8.3\times 10^{3}$
M16 | $M_{\mathrm{star,3pc}}=3\times 10^{6}{\rm M}_{\odot}$ | $0.01$ | $2\times 10^{-4}$ | $6\times 10^{-4}$ | 0.98 | $0.01$ | $5\times 10^{-4}$ | 0 | $5.3\times 10^{2}$
M17 | $\gamma_{\mathrm{\rho,BH}}=1.5$ | $1\times 10^{-3}$ | $0$ | $7\times 10^{-5}$ | 0.994 | $4\times 10^{-3}$ | $2\times 10^{-4}$ | 0 | $9.8\times 10^{2}$
M18 | ${\beta}_{\mathrm{v,BH}}=0.4$ | $9\times 10^{-3}$ | $1\times 10^{-3}$ | $2\times 10^{-3}$ | 0.97 | $0.01$ | $1\times 10^{-3}$ | 0 | $1.2\times 10^{3}$
M19 | ${\beta}_{\mathrm{v,NS}}=0.2$ | $6\times 10^{-3}$ | $7\times 10^{-4}$ | $8\times 10^{-4}$ | 0.987 | $6\times 10^{-3}$ | $7\times 10^{-4}$ | 0.06 | $2.2\times 10^{3}$
M20 | BS in 2D | $3\times 10^{-3}$ | $1\times 10^{-4}$ | $6\times 10^{-4}$ | 0.987 | $9\times 10^{-3}$ | $3\times 10^{-4}$ | 0.07 | $4.6\times 10^{3}$
M21 | BS in $2\rightarrow 3$D | $3\times 10^{-3}$ | $2\times 10^{-4}$ | $4\times 10^{-4}$ | 0.991 | $6\times 10^{-3}$ | $3\times 10^{-4}$ | 0.05 | $3.4\times 10^{3}$
M22 | $m_{\mathrm{1gBH}}=10\,{\rm M}_{\odot}$ | $3\times 10^{-3}$ | $3\times 10^{-4}$ | $3\times 10^{-4}$ | 0.991 | $5\times 10^{-3}$ | $4\times 10^{-4}$ | $0.01$ | $1.6\times 10^{3}$
M1 | $t=1\,\mathrm{Myr}$ | $4\times 10^{-3}$ | $3\times 10^{-4}$ | $9\times 10^{-4}$ | 0.986 | $8\times 10^{-3}$ | $8\times 10^{-4}$ | $0$ | $2.6\times 10^{2}$
M1 | $t=100\,\mathrm{Myr}$ | $0.02$ | $9\times 10^{-4}$ | $1\times 10^{-3}$ | 0.96 | $0.02$ | $8\times 10^{-4}$ | $8\times 10^{-3}$ | $7.3\times 10^{3}$
M12 | ${\dot{M}}_{\mathrm{out}}=0.01{\dot{M}}_{\mathrm{Edd}},t=100\,\mathrm{Myr}$ | $7\times 10^{-3}$ | $5\times 10^{-4}$ | $1\times 10^{-3}$ | 0.98 | $7\times 10^{-3}$ | $7\times 10^{-4}$ | 0.05 | $4.6\times 10^{3}$
M12 | ${\dot{M}}_{\mathrm{out}}=0.01{\dot{M}}_{\mathrm{Edd}},t=1\,\mathrm{Gyr}$ | $4\times 10^{-3}$ | $2\times 10^{-3}$ | $6\times 10^{-3}$ | 0.92 | $0.03$ | $3\times 10^{-3}$ | 0.07 | $1.7\times 10^{4}$
M13 | ${\dot{M}}_{\mathrm{out}}=0.001{\dot{M}}_{\mathrm{Edd}},t=100\,\mathrm{Myr}$ | $3\times 10^{-3}$ | $4\times 10^{-4}$ | $5\times 10^{-4}$ | 0.993 | $3\times 10^{-3}$ | $4\times 10^{-4}$ | $4\times 10^{-3}$ | $2.2\times 10^{3}$

##  Appendix A Alternative models

In this section, we present the parameter dependence of the results using a
total of 22 different models, listed in Table 2. We also discuss merging mass
distributions in the AGN channel.

##  Appendix B Mergers in upper-mass gap

We find that the fraction of mergers of massive binaries of $\geq 130\,{\rm
M}_{\odot}$ ($f_{\mathrm{M130}}$) is zero when the BH number in inner regions
($r\lesssim 0.01\,\mathrm{pc}$) of the AGN disk is low due to the small size
of the AGN disk (model M4), the high velocity dispersion of BHs (models M5,
M16), no radial migration (model M7), low mass for NSCs (model M14), shallow
density profile of BHs (model M15), or short duration of an AGN phase
($t=1\,\mathrm{Myr}$, model M1). On the other hand, $f_{\mathrm{M130}}$ is
enhanced if the initial masses of 1g BHs are massive (model M2), super-
Eddington accretion is allowed (model M3), large number of BHs initially
exists (model M13), or BHs can efficiently migrate to inner regions of
$r\lesssim 10^{-3}\,\mathrm{pc}$. In models M8–M11, the fraction of mergers in
the inner regions ($r\lesssim 10^{-3}\,\mathrm{pc}$) is higher than that in
model M1. In the inner regions, massive BHs can merge with high binary
eccentricity due to GW-capture binary formation (Paper III, see also Gondán &
Kocsis 2020) which may be consistent with possible high eccentricity in
GW190521 implied by Gayathri et al. (2020).

Here, the maximum masses of merging binaries in the fiducial model ($\sim
185\,{\rm M}_{\odot}$) are significantly reduced from those in Papers I ($\sim
10^{4}\,{\rm M}_{\odot}$). This is due to the exchange interactions and
disruption of softer binaries during binary-binary interactions, which are not
included previously. With the exchanges included, both the number and maximum
mass of mergers become smaller, because massive objects are retained in
binaries after the exchanges and the more energy must be extracted for the
more massive binaries to merge, during which many lighter COs receive strong
kicks, reducing the number of COs in AGN disks.

Also, in our fiducial model, mergers between two high-generation BHs are
common and the majority of them have high mass-ratio close to unity, as is the
case for GW190521. Such high-g-high-g mergers do not occur in the model by
Yang et al. (2019), who assumed that mergers occurs in the migration traps
shortly after being captured in the disk. In that case, runaway growth of a
single object by $n$g-$1$g mergers took place instead. In AGN disks outside
the migration traps, high-g-high-g mergers can be more frequent because of the
following reasons: (i) it takes longer for binary hardening and merger, (ii)
objects kicked out from the disk by BS interactions interact less with others,
and (iii) different objects orbit and merge at different distances from the
SMBH due to variations in initial positions, gap-opening regions, and
migration speeds. Combinations of these enable objects to avoid merging into
one massive object in contrary to mergers in the migration traps.

In the fiducial model, at $t=1$, $10$, and $100~{}\mathrm{Myr}$, respectively,
the maximum masses of merging binaries are $130$, $185$, and $185~{}{\rm
M}_{\odot}$, and the numbers of mergers are $260$, $1,600$, and $7,300$, if
the AGN disk exists for these durations. Thus, the merger rate decreases with
time, as the rate of new BH captures in the disk by gas dynamical friction
decreases. This is because of the strong dependence of the capture timescale
on the initial CO velocity relative to the local AGN disk gas motion ($\propto
v_{\mathrm{ini}}^{4}$, Eq. 26 in paper I). COs having low $v_{\mathrm{ini}}$
are quickly captured, while those with high $v_{\mathrm{ini}}$ are typically
not captured by the AGN disk within $100\,\mathrm{Myr}$. Also, the maximum
masses are not sensitive to disk lifetimes beyond $t\geq 10\,\mathrm{Myr}$.
This is because mergers become infrequent and massive BHs often migrate to the
inner boundary. Here, the number of BHs that migrate to the inner boundary of
$r=10^{-4}\,\mathrm{pc}$ is $\sim 90$ at $t=100\,\mathrm{Myr}$, while it is
$\sim 10$ at $t=10\,\mathrm{Myr}$. Note that objects which reach $r\leq
10^{-4}\,\mathrm{pc}$ are not followed in our model, and are instead removed
from the simulation. In these regions, COs may merge with each other if they
accumulate in migration traps (Bellovary et al., 2016) or they may accrete
onto the SMBH and observed by LISA as the migration timescale due to GW
emission is comparable to the AGN lifetime ($\sim
30\,\mathrm{Myr}(r/10^{-4}\,\mathrm{pc})^{4}(M_{\mathrm{SMBH}}/4\times
10^{6}\,{\rm M}_{\odot})^{-2}(m_{\mathrm{CO}}/100\,{\rm M}_{\odot})^{-1}$,
where $M_{\mathrm{SMBH}}$ and $m_{\mathrm{CO}}$ are the SMBH and the CO
masses). Note that as massive BHs are rare, the evolution of maximum masses
likely suffers some variations. For example, with an independent realization
of the initial condition, the maximum masses of merging binaries are $70$,
$138$, and $186~{}{\rm M}_{\odot}$, at $t=1$, $10$, and $100~{}\mathrm{Myr}$,
respectively. In all cases, the growth of the maximum mass is slowed down in
later phases.

In the fiducial model, $f_{\mathrm{M130}}$ decreases in $t\geq
10\,\mathrm{Myr}$ due to inefficient capture of BHs in later phases and
migration of massive BHs to the inner boundary of the AGN disk as discussed
above. On the other hand, $f_{\mathrm{M130}}$ increases in $t\geq
10\,\mathrm{Myr}$ in model M12 and M13, in which the gas accretion rate at the
outer boundary is low (${\dot{M}}_{\mathrm{out}}=0.01$ and
$0.001\,{\dot{M}}_{\mathrm{Edd}}$). In addition, in model M12, the maximum
mass of merging binaries increases with time as those are $46$, $142$, $228$,
and $283~{}{\rm M}_{\odot}$ at $t=1$, $10$, $100$, and $1000~{}\mathrm{Myr}$,
respectively. Thus, the maximum merging mass and $f_{\mathrm{M130}}$ evolve
more in later phases for low ${\dot{M}}_{\mathrm{Edd}}$ models compared to
that in the fiducial model. These are due to the difference of the gas density
of the AGN disk. Inefficient gaseous processes (capture to the AGN disk and
migration) due to low gas density facilitate growth of BHs in later phases.
Such low ${\dot{M}}_{\mathrm{out}}$ cases may be more important for
understanding properties of mergers in the AGN channel as the number of
mergers from AGNs with low Eddington rates may dominate over those with high
Eddington rates depending on the Eddington rate distribution of AGNs (Hopkins
et al., 2005; Greene & Ho, 2007; Shankar et al., 2013), which is uncertain. If
the merger rate is dominated by AGNs with low Eddington rates, the rate is
significantly influenced by the transition radius from a standard disk to
advection-dominated accretion flow, which possibly comes to $\gtrsim
10^{-3}$–$10^{-2}\,\mathrm{pc}$ and affect merger processes for
${\dot{M}}_{\mathrm{out}}\lesssim 10^{-3}\,{\dot{M}}_{\mathrm{Edd}}$ (Eq. 3 of
Czerny et al. 2019), although the typical value of the radius is uncertain
(e.g. Czerny et al., 2019).

##  Appendix C Mergers in lower mass gap

The fraction of BH-LGO mergers is enhanced by several effects. In model M4, we
assume that COs initially exist within
$r_{\mathrm{out,BH}}=0.03\,\mathrm{pc}$. This resembles the situation in which
the size of the AGN disk is $0.03\,\mathrm{pc}$. In this model,
$f_{\mathrm{BHLGO}}$ is enhanced to $0.10$. In the inner regions of the AGN
disk, the capture timescale of COs by the AGN disk is shorter. As a large
fraction of NSs and BHs are rapidly captured, the ratio of the number of NSs
over BHs in the AGN disk is high and accordingly the rates of NS-NS and BH-LGO
mergers are also relatively high.

In model M5, in which the initial velocity dispersion of BHs and NSs is virial
($\beta_{\mathrm{v,BH}}=\beta_{\mathrm{v,NS}}=1$), $f_{\mathrm{BHLGO}}$ is
enhanced to $0.03$. As the initial velocity dispersion of BHs and NSs
increases, the numbers of BHs and NSs in the AGN disks decrease. Then, binary
NSs have a high probability to merge without interacting with BHs and
experiencing exchange of binary components. The high velocity dispersion of
BHs might be realized if the AGN phase begins before dynamics of NSCs is
regulated by vector resonant relaxation e.g. due to randomization of orbits of
COs by major merger of galaxies, or if NSCs have negligible rotation on
average.

In models M6 and M7, we assume that NSs have already merged prior to the AGN
phase. Such situation may be realized if the high rate of NS-NS mergers ($\sim
1000\,\mathrm{Gpc^{-3}yr^{-1}}$) inferred from GW observations are driven by
isolated binary evolution (but see Belczynski et al., 2018) and such mergers
are not prevented by soft-BS interactions in NSCs, or NS mergers are
facilitated by the Kozai-Lidov mechanism (e.g. Fragione et al., 2019). To
investigate the effects of this possibility, we assume that a fraction of NSs
($f_{\mathrm{LGO}}$) already merged with NSs prior to the start of the AGN
episode, and LGOs may have initial masses with $2.6\,{\rm M}_{\odot}$. We
investigate a model with $f_{\mathrm{LGO}}=0.1$ (M6), which would present an
upper limit as the rate of core-collapse supernovae is $\sim
10^{5}\,\mathrm{Gpc^{-3}yr^{-1}}$ (Madau & Dickinson, 2014) and that of NS-NS
mergers is $\sim 100$–$4000\,\mathrm{Gpc^{-3}yr^{-1}}$ (Abbott et al., 2019).
$f_{\mathrm{BHLGO}}$ is moderately enhanced to $0.01$ for
$f_{\mathrm{LGO}}=0.1$ from its value of $4\times 10^{-3}$ for
$f_{\mathrm{LGO}}=0$. Thus, isolated binary evolution may contribute to the
fraction of mergers in the lower mass gap ($f_{\mathrm{BHLGO}}$) in the AGN
channel.

In model M8, in which COs do not migrate toward the SMBH, $f_{\mathrm{BHLGO}}$
is higher as the probability that the NS binaries interact with BHs and the
binary components are exchanged is reduced. The efficiency of migration of COs
in AGN disks is uncertain, while it affects various properties of mergers in
AGN disks. Thus, this process is worth investigating in more detail in future
work to better understand the implications for mergers in AGN disks.

Figure 4.—  The distribution of the chirp (left) and primary (right) masses
for model M1 at $t=10\,\mathrm{Myr}$ (panels a and b) and model M22 at
$t=10\,\mathrm{Myr}$ (panels c and d) and $t=30\,\mathrm{Myr}$ (panels e and
f).

##  Appendix D Merging mass distribution

Here, we discuss the merging mass distribution. As a large number of events
are reported in LIGO/Virgo O3a (Abbott et al., 2020c), we can assess the
consistency between predicted and observed distributions more precisely.
Recently, The LIGO Scientific Collaboration et al. (2020c) and Tiwari &
Fairhurst (2020) analyzed the LIGO/Virgo’s O1-O3a data and reconstructed
merging mass distributions. The LIGO Scientific Collaboration et al. (2020c)
reconstructed the primary mass distribution by a superposition of a (broken)
power law and $0$–$2$ Gaussian distributions, while Tiwari & Fairhurst (2020)
reconstructed the chirp mass distribution by a superposition of twelve
Gaussian distributions, which can flexibly recover features of the mass
distribution. The LIGO Scientific Collaboration et al. (2020c) suggested that
the $m_{1}$ distribution is more consistent with a broken power law with a
break at $m_{1}\sim 40\,{\rm M}_{\odot}$, or a power law with a Gaussian
feature peaking at $m_{1}\sim 34\,{\rm M}_{\odot}$. Interestingly, Tiwari &
Fairhurst (2020) identified four peaks in the chirp mass distribution at 8,
14, 26, and $45\,{\rm M}_{\odot}$, which are roughly separated by a factor of
around two, possibly suggesting frequent hierarchical mergers. Here, high
values for $\chi_{\mathrm{p}}$ inferred for observed GW events (Gerosa et al.,
2020a) may also support this scenario. As hierarchical mergers are presumably
most frequent in the AGN channel, it would be interesting to check whether
such features can be reproduced in our model.

In model M1, peak-like structures are less clear in the chirp mass
distribution (Fig. 4 a) despite of frequent hierarchical mergers. On the other
hand, as the merging mass distribution is significantly affected by the
adopted initial mass function of BHs, we performed model M22 in which the
initial masses of all 1g BHs are set to $10\,{\rm M}_{\odot}$. Such a spiky
initial mass function may be allowed depending on prescriptions of stellar
mass loss (Belczynski et al., 2010). In this model, peaks in the
$m_{\mathrm{chirp}}$ distribution become much clearer (panels c and e),
although the number of peaks is larger than that implied by Tiwari & Fairhurst
(2020), due to the contribution of mergers among various generation BHs like
1g-2g or 2g-3g. This model predicts that further peaks might be detectable by
observing a larger number of events. On the other hand, the contribution of
mergers among 1g-2g or 2g-3g BHs is supposedly influenced by frequency of BS
interactions, which depends on several parameters and processes, such as
$\delta_{\mathrm{IMF}}$, $\beta_{\mathrm{v,BH}}$, and the efficiency of
migration due to gaseous torque. Also, some peaks become smooth and close to a
power-law like shape if we adopt different initial mass functions for BHs as
in model M1 (panel a) or as time passes (panels c and e). Observational errors
on the chirp mass further make the merging mass distribution somewhat
smoother. As properties of peaks are influenced by several parameters,
uncertainties, and observational effects, the features inferred by Tiwari &
Fairhurst (2020) might be reproduced.

Meanwhile, the $m_{1}$ distributions (Fig. 4 b, d, and e) look roughly similar
to those reconstructed by The LIGO Scientific Collaboration et al. (2020c).
Note that the errors on $m_{1}$ are typically much larger than that on the
chirp mass, which presumably reduces peaks in the observed $m_{1}$
distribution (Fig. 4 d and f). Hence, there is a possibility that the AGN
channel can explain the observed features of merging masses. To constrain the
contribution from the AGN channel, quantitative comparisons would be required,
which will be conducted in a following work.

Table 3Fiducial values of our model parameters.  Parameter | Fiducial value
---|---
Spacial directions in which BS interactions occur | isotropic in 3D
Number of temporary binary COs formed during a BS interaction | $N_{\mathrm{int}}=20$
Initial BH spin magnitude | $|{\bm{a}}|=0$
Angular momentum directions of circum-CO disks |  $\hat{{\bm{J}}}_{\mathrm{CCOD}}=\hat{{\bm{J}}}_{\mathrm{AGN}}$ for single COs,
|  $\hat{{\bm{J}}}_{\mathrm{CCOD}}=\hat{{\bm{J}}}_{\mathrm{bin}}$ for COs in
binaries
Ratio of viscosity responsible for warp propagation |
over that for transferring angular momentum | $\nu_{2}/\nu_{1}=10$
Alignment efficiency of the binary orbital angular momentum due to gas capture | $f_{\mathrm{rot}}=1$
(Eq. 14 in Paper II) |
Mass of the central SMBH | $M_{\mathrm{SMBH}}=4\times 10^{6}\,{\rm M}_{\odot}$
Gas accretion rate at the outer radius of the simulation ($5\,\mathrm{pc}$) |  ${\dot{M}}_{\mathrm{out}}=0.1\,{\dot{M}}_{\mathrm{Edd}}$ with $\eta=0.1$
Fraction of pre-existing binaries | $f_{\mathrm{pre}}=0.15$
Power-law exponent for the initial density profile for BHs and NSs |  $\gamma_{\mathrm{\rho,BH}}=0$, $\gamma_{\mathrm{\rho,NS}}=-0.5$
Initial velocity anisotropy parameter |
such that $\beta_{\mathrm{v,BH}}v_{\mathrm{kep}}(r)$ is the BH velocity dispersion |  $\beta_{\mathrm{v,BH}}=0.2$, $\beta_{\mathrm{v,NS}}=0.4$
Efficiency of angular momentum transport in the $\alpha$-disk | $\alpha_{\mathrm{SS}}=0.1$
Stellar mass within 3 pc | $M_{\mathrm{star,3pc}}=10^{7}\,{\rm M}_{\odot}$
Stellar initial mass function slope | $\delta_{\mathrm{IMF}}=2.35$
Angular momentum transfer parameter in the outer star forming regions | $m_{\mathrm{AM}}=0.15$
(Eq. C8 in Thompson et al. 2005) |
Accretion rate in Eddington units onto |
stellar-mass COs with a radiative efficiency $\eta=0.1$ | $\Gamma_{\mathrm{Edd,cir}}=1$
Numerical time-step parameter | $\eta_{t}=0.1$
Number of radial cells storing physical quantities | $N_{\mathrm{cell}}=120$
Maximum and minimum $r$ for the initial CO distribution |  $r_{\mathrm{in,BH}}=10^{-4}$ pc, $r_{\mathrm{out,BH}}=3$ pc
Initial number of BHs and NSs within 3 pc |  $N_{\mathrm{BH,ini}}=2\times 10^{4}$, $N_{\mathrm{NS,ini}}=5\times 10^{4}$


redacted 001

# Model Predictive Simulation Using Structured Graphical Models and
Transformers

Xinghua Lou  Google DeepMind  Meet Dave  Google DeepMind  Shrinu Kushagra
Google DeepMind  Miguel Lázaro-Gredilla  Google DeepMind  Kevin Murphy  Google
DeepMind

###### Abstract

We propose an approach to simulating trajectories of multiple interacting
agents (road users) based on transformers and probabilistic graphical models
(PGMs), and apply it to the Waymo SimAgents challenge. The transformer
baseline is based on the MTR model (Shi et al., 2024), which predicts multiple
future trajectories conditioned on the past trajectories and static road
layout features. We then improve upon these generated trajectories using a
PGM, which contains factors which encode prior knowledge, such as a preference
for smooth trajectories, and avoidance of collisions with static obstacles and
other moving agents. We perform (approximate) MAP inference in this PGM using
the Gauss-Newton method. Finally we sample $K=32$ trajectories for each of the
$N\sim 100$ agents for the next $T=8\Delta$ time steps, where $\Delta=10$ is
the sampling rate per second. Following the Model Predictive Control (MPC)
paradigm, we only return the first element of our forecasted trajectories at
each step, and then we replan, so that the simulation can constantly adapt to
its changing environment. We therefore call our approach "Model Predictive
Simulation" or MPS. We show that MPS improves upon the MTR baseline,
especially in safety critical metrics such as collision rate. Furthermore, our
approach is compatible with any underlying forecasting model, and does not
require extra training, so we believe it is a valuable contribution to the
community.

###### keywords:

Multi-Agent Planning, Probabilistic Graphical Model, Model Predictive Control,
Transformer

##  1 Introduction

The use of transformers to create generative models to simulate agent
trajectories, trained on large datasets such as Waymo Open Data (Ettinger et
al., 2021), has become very popular in recent years. Most previous work has
been focusing on improving the architecture (Nayakanti et al., 2023; Shi et
al., 2024), the training objective (Ngiam et al., 2021; Shi et al., 2024), the
trajectory representation (Seff et al., 2023; Philion et al., 2023) or the
speed (Zhou et al., 2023) of these transformer-based models.

This paper tackles the problem from an orthogonal and complementary angle –
namely the use of prior knowledge, encoded using a probabilistic graphical
model (PGM). We perform approximate MAP inference in the PGM to “post process”
the trajectory proposals from a base transformer model, to increase their
realism and compliance with constraints, such as collision avoidance.

To ensure that our predicted forecasts are adaptive to the changing
environment, we replan at each step, following the principle of model
predictive control (MPC), which is widely used for controlling complex
dynamical systems (Schwenzer et al., 2021). We therefore call our approach
Model Predictive Simulation (MPS).

Our MPS approach differs from previous PGM methods for trajectory simulation,
such as JFP (Luo et al., 2023), in several ways. First, we explicitly include
(data-dependent) factors for collision avoidance and smooth trajectories, so
we have better control over the generated trajectories. Second, our approach
is iterative (being based on MPC), while JFP commits to the trajectory
proposals at $t=0$ and is thus open loop. Third, our approach uses the Gauss-
Newton method to compute the joint MAP estimate, whereas JFP is based on
discrete belief propagation methods to choose amongst a finite set of
candidate trajectories.

##  2 Method

### Outer loop

Input: Scene context $\mathbf{c}$, num. agents $N$, num. samples $K$,
trajetory length $T$

Output: Sampled trajectories, $h_{1:N}^{1:K,1:T}$

for  _$k=1$ to $K$_ do

$s_{1:N}^{k,0}=\text{init-trajectory}(\mathbf{c})$

for  _$t=1$ to $T$_ do

Sample $r_{1:N}^{k,t}=\text{MPS}(\mathbf{c},s_{1:N}^{k,1:t-1})$

Extend $s_{1:N}^{k,1:t}=\text{append}(s_{1:N}^{k,1:t},r_{1:N}^{k,t})$

end for

end for

Algorithm 1 SimAgents outer loop

The overall simulation pseudocode is shown in Algo. 1. It generates a set of
$K=32$ trajectories, each of length $T=80$, for $N$ agents given the scene
context $\mathbf{c}$. (The exact value of $N$ depends on the number of agents
that are visible in $\mathbf{c}$.) We denote the generated output by
$s_{1:N}^{1:K,1:T}$, where $s_{i}^{k,t}=[x,y,\dot{x},\dot{y}]$ is the state
(2d location and velocity) of the $i$’th agent in sample $k$.

### Inner loop

Input: Scene context $\mathbf{c}$, agent history $h_{1:N}$, num. agents $N$,
future planning horizon $F$, number of rollouts $J$, transformer proposal
$\pi$

Output: Predicted next state for each agent, $r_{1:N}$

for  _$j=1$ to $J$_ do

Sample $(a_{1:N}^{j,1:F},g_{1:N}^{j})\sim\pi(\mathbf{c},h_{1:N},F)$

$G^{j}=\text{BuildFactorGraph}(a_{1:N}^{j,1:F},g_{1:N}^{j},\mathbf{c})$

Initialize $s_{1:N}^{j,1:F}=a_{1:N}^{j,1:F}$

$(s_{1:N}^{j,1:F},E^{j})=\text{Inference}(G^{j},s_{1:N}^{j,1:F})$

end for

Sample $j^{*}\sim\text{SoftMin}(E^{1:J})$

return $s_{1:N}^{j^{*},1}$

Algorithm 2 Model Predictive Simulation

At each step $t$, the simulator calls our MPS algorithm to generate a
prediction for the next state of each agent. The pseudocode for this is shown
in Algo. 2. The approach is as follows. First we use the MTR transformer model
$\pi$ (Shi et al., 2024) to sample a set of $N$ goal locations, $g_{1:N}^{j}$,
one for each agent, as well as a sequence of anchor points leading to each
goal, $a_{1:N}^{j,1:F}$, where $F$ is the planning or forecast horizon. We do
this $J=60$ times in parallel, to create a set of possible futures. We then
use the PGM to generate $J$ joint trajectories (for all $N$ agents), using the
method described below. Finally we evaluate the energy of each generated
trajectory, $E^{j}$, sample one of the low energy (high probability) ones to
get $s_{1:N}^{j^{*},1:F}$, and return the first step of this sampled
trajectory, $s_{1:N}^{j^{*},1}$.

### Graphical model

Figure 1: Factor Graph for $N=2$ agents unrolled for $T$ planning steps.
Circles are random variables, gray squares are fixed factors.  $\displaystyle
p(\mathbf{s}^{1:F}_{1..N}\mid\mathbf{c},\mathbf{a}_{1:N}^{1:F},g_{1:N})\propto$
$\displaystyle\prod_{i=1}^{N}\left[f_{G}(s_{i}^{F}|g_{i})\cdot\prod_{t=1}^{F-1}f_{M}(s_{i}^{t}|a_{i}^{t})\right]\cdot$
$\displaystyle\prod_{i=1}^{N}\left[\prod_{t=2}^{F}f_{L}(s_{i}^{t-1},s_{i}^{t})\cdot
f_{A}(s_{i}^{t-1},s_{i}^{t})\right]\cdot$
$\displaystyle\prod_{i=1}^{N}\prod_{t=1}^{F}f_{O}(s_{i}^{t}\mid\mathbf{c})\cdot\prod_{i\neq
j}\prod_{t=1}^{F}f_{C}(s_{i}^{t},s_{j}^{t})$ Figure 2: Joint probability
model.

The key to our method is the probabilistic graphical model (PGM) for improving
upon the proposed trajectories by MTR. The factor graph is shown in Fig. 1 and
the corresponding conditional joint distribution is given in Fig. 2. The model
was inspired by (Patwardhan et al., 2022) who uses Gaussian belief
propagation. We now explain each of the factors.

First we have factors which compare a candidate trajectory to the original
proposal. The motion factor is defined as
$f_{\mathrm{M}}(\mathbf{s}_{i}^{t})=|s_{i}^{t}-a_{i}^{t}|$, where $a_{i}^{t}$
is the predicted location (anchor point) for agent $i$ at time $t$ as computed
by $\pi$. This ensures the trajectory stays close to the initial proposal. The
proximity to goal factor is defined as
$f_{\mathrm{G}}(s_{i}^{F})=|s_{i}^{F}-g_{i}|$, where $g_{i}$ is the goal for
agent $i$ predicted by $\pi$. This ensures the trajectory ends close to where
we expect.

Second we have factors defined from "physics". We define a factor that
penalizes deviation from linear motion:
$f_{\mathrm{L}}(s,s^{\prime})=|s^{\prime}_{xy}-(s_{xy}+s_{\dot{x}\dot{y}}\Delta_{t})|$,
where $s_{xy}$ are the location components of $s$, $s_{\dot{x}\dot{y}}$ are
the velocity components of $s$, and $\Delta t$ is the sampling rate. We also
define a factor that penalizes change in direction:
$f_{\mathrm{A}}(s,s^{\prime})=|s_{\dot{x}\dot{y}}-s^{\prime}_{\dot{x}\dot{y}}|$.
We used weight 2.0 for $f_{\mathrm{A}}$ and 1.0 for all other factors.

Third we have factors derived from static obstacles on the road:
$f_{\mathrm{O}}(s\mid\mathbf{c})=\max_{(x,y)\in\mathbf{c}_{\mathrm{RE}}}\mathrm{G}(x,y\mid
s)$, where $\mathbf{c}_{\mathrm{RE}}$ represents the coordinates of the road
edges (part of the context $\mathbf{c}$) and $\mathrm{G}(x,y\mid s)$ is a
Gaussian field centered and rotated according to the agent’s location
$s_{xy}$.

Finally, we have pairwise collision factors between agents:
$f_{\mathrm{C}}(s,s^{\prime})=\max_{(x,y)\in\mathrm{CCP}(s^{\prime})}\mathrm{G}(x,y\mid
s)$, where $\mathrm{G}(x,y\mid s)$ is a Gaussian field for agent $s$, and
$\mathrm{CCP}(s^{\prime})$ are the 9 collision checking points (CCP) for the
other agent $s^{\prime}$ (4 corners, 4 centers of the sides, and center of the
agent).

### Inference

Inference on the factor graph is equivalent to minimizing a non-linear, non-
convex quadratic optimization problem defined over $\mathbf{s}_{1:N}^{1:F}$.
For efficiency reasons, we developed a two-step approach. First, we use the
Gauss–Newton method to solve a partial model that only consists of
$f_{\mathrm{M}}$, $f_{\mathrm{L}}$ and $f_{\mathrm{A}}$ factors, as these can
all be evaluated in parallel across agents using $N$ individual trajectory
models. This step produces smoothed trajectories, which are then frozen.
Second, we sample joint trajectories for agents according to their probability
(unnormalized energy), and use the $f_{\mathrm{O}}$ and $f_{\mathrm{C}}$
factors to score their quality. After repeating this $J$ times, the best joint
trajectories are sampled from a softmin operation over the scores of the $J$
samples.

WAYMO | META METRIC | KINEMATIC | INTERACTIVE | MAP |
---|---|---|---|---|---
LEADERBOARD | REALISM |  | LINEAR
---
SPEED
| LINEAR
---
ACCEL.
| ANG.
---
SPEED
| ANG.
---
ACCEL.
| DIST.
---
TO OBJ.
COLLISION | TTC |  | DIST.
---
TO ROAD
OFFROAD | minADE $\downarrow$
SMART | 0.7511 | 0.3646 | 0.4057 | 0.4231 | 0.5844 | 0.3769 | 0.9655 | 0.8317 | 0.6590 | 0.9362 | 1.5447
MVTE | 0.7301 | 0.3506 | 0.3530 | 0.4974 | 0.5999 | 0.3742 | 0.9049 | 0.8309 | 0.6655 | 0.9071 | 1.6769
MPS (Ours) | 0.7416 | 0.3137 | 0.3049 | 0.4705 | 0.5834 | 0.3593 | 0.9629 | 0.8070 | 0.6651 | 0.9366 | 1.4841

Table 1: WOSAC Leaderboard: SMART (2024 winner) Vs. MVTE (2023 winner) Vs. MPS
(ours).

##  3 Experimental Evaluation

Benchmark We evaluated MPS on the 2024 Waymo Open Sim Agents Challenge
(Montali et al., 2023), where the task is simulating 32 realistic rollouts of
all agents in the scene given their 1s history for 8s into the future. The
simulation needs to be closed-loop and factorized between the ADV and other
agents, which MPS satisfies naturally.

Implementation Details We implemented the factors and the inference in JAX
111https://github.com/google/jax and JAXopt 222https://jaxopt.github.io/ for
the Gauss-Newton method. We leveraged JAX’s just-in-time (JIT) compilation and
observed great scalability. For speed up, we take 10 immediate next steps at
each MPS iteration.

We trained our own MTR model $\pi$ using the open source code
333https://github.com/sshaoshuai/MTR . We removed local attention and reduced
the source polylines to 512. The training data is augmented by adding extra
interacting agents, and by applying random history dropouts. We followed the
original training setup except the number of epochs (50), the batch size (8)
and the LR schedule ([25, 30, 35, 40, 45]). Training took about 3 days on 16
A100s.

We used only the official Waymo Open Motion Dataset v1.2.1 and did not use any
Lidar or Camera data. We did not need any additional training and we did not
use ensembles.

Sim Agents 2024 Results We ranked number 4 among all methods (Table 1). We
outperformed the 2023 winner MVTE (Wang et al., 2023) which also uses MTR (Shi
et al., 2024), and are approximately 1 point behind the 2024 winner SMART (Wu
et al., 2024). MPS achieved near-top performance in a few safety critical
metrics such as COLLISION and OFFROAD, showing the effectiveness of the priors
in our model. MPS showed a lack of performance in LINEAR SPEED / ACCEL. We
speculate this is because MPS can generate diverse rollouts that are very
different from the logged data used for metric evaluation.

| META METRIC | KINEMATIC | INTERACTIVE | MAP |
---|---|---|---|---|---
METHOD | REALISM |  | LINEAR
---
SPEED
| LINEAR
---
ACCEL.
| ANG.
---
SPEED
| ANG.
---
ACCEL.
| DIST.
---
TO OBJ.
COLLISION | TTC |  | DIST.
---
TO ROAD
OFFROAD | minADE $\downarrow$
MTR+RAND | 0.7019 | 0.3922 | 0.3530 | 0.3899 | 0.3304 | 0.3691 | 0.8491 | 0.8164 | 0.6706 | 0.9207 | 1.3084
MPS | 0.7418 | 0.3158 | 0.3056 | 0.4664 | 0.5818 | 0.3604 | 0.9617 | 0.8094 | 0.6651 | 0.9374 | 1.4841

Table 2: Abalation study – comparing MPS to MTR with random trajectory
sampling.

Ablation Study To evaluate the value of the PGM priors, we compare MPS to the
same MTR model with random trajectory sampling (MTR+RAND) on the validation
dataset. As shown in Table 2, MPS improved safety-critical metrics such as
COLLISION, OFFROAD and the overall REALISM score, while lacked performance at
LINEAR SPEED / ACCEL for the same reason discussed above.

Qualitative Study Qualitatively, MPS generates diverse (multi-modal)
predictions (Fig. 3), and each prediction contains realistic traffic patterns
such as lane merging, unprotected left turn, yielding, among others (Fig. 4).

Figure 3: MPS creates diverse (multi-modal) rollouts. Figure 4: MPS creates
realistic traffic patterns.

##  4 Conclusion

We explored an approach that can improve on any trajectory simulation model by
adding domain-specific priors, and performing inference in the corresponding
PGM. We believe combing prior-driven (top-down) and data-driven (bottom-up)
methods is key to building robust and reliable autonomous driving planning and
simulation.444We thank Joseph Ortiz and Wolfgang Lehrach for many useful
discussions and suggestions.

*

## References

  * Ettinger et al. (2021) S. Ettinger, S. Cheng, B. Caine, C. Liu, H. Zhao, S. Pradhan, Y. Chai, B. Sapp, C. R. Qi, Y. Zhou, et al.  Large scale interactive motion forecasting for autonomous driving: The waymo open motion dataset.  In _ICCV_ , 2021.
  * Luo et al. (2023) W. Luo, C. Park, A. Cornman, B. Sapp, and D. Anguelov.  Jfp: Joint future prediction with interactive multi-agent modeling for autonomous driving.  In _CoRL_ , 2023.
  * Montali et al. (2023) N. Montali, J. Lambert, P. Mougin, A. Kuefler, N. Rhinehart, M. Li, C. Gulino, T. Emrich, Z. Yang, S. Whiteson, B. White, and D. Anguelov.  The waymo open sim agents challenge.  In _NIPS Datasets Track_ , May 2023.
  * Nayakanti et al. (2023) N. Nayakanti, R. Al-Rfou, A. Zhou, K. Goel, K. S. Refaat, and B. Sapp.  Wayformer: Motion forecasting via simple & efficient attention networks.  In _ICRA_. IEEE, 2023.
  * Ngiam et al. (2021) J. Ngiam, B. Caine, V. Vasudevan, Z. Zhang, H.-T. L. Chiang, J. Ling, R. Roelofs, A. Bewley, C. Liu, A. Venugopal, et al.  Scene transformer: A unified architecture for predicting multiple agent trajectories.  _arXiv preprint arXiv:2106.08417_ , 2021.
  * Patwardhan et al. (2022) A. Patwardhan, R. Murai, and A. J. Davison.  Distributing collaborative Multi-Robot planning with gaussian belief propagation.  In _IEEE Robotics and Automation Letters_ , Mar. 2022.
  * Philion et al. (2023) J. Philion, X. B. Peng, and S. Fidler.  Trajeglish: Learning the language of driving scenarios.  _arXiv preprint arXiv:2312.04535_ , 2023.
  * Schwenzer et al. (2021) M. Schwenzer, M. Ay, T. Bergs, and D. Abel.  Review on model predictive control: An engineering perspective.  _The International Journal of Advanced Manufacturing Technology_ , 117(5):1327–1349, 2021.
  * Seff et al. (2023) A. Seff, B. Cera, D. Chen, M. Ng, A. Zhou, N. Nayakanti, K. S. Refaat, R. Al-Rfou, and B. Sapp.  Motionlm: Multi-agent motion forecasting as language modeling.  In _ICCV_ , 2023.
  * Shi et al. (2024) S. Shi, L. Jiang, D. Dai, and B. Schiele.  Mtr++: Multi-agent motion prediction with symmetric scene modeling and guided intention querying.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2024.
  * Wang et al. (2023) Y. Wang, T. Zhao, and F. Yi.  Multiverse transformer: 1st place solution for waymo open sim agents challenge 2023.  _arXiv preprint arXiv:2306.11868_ , 2023.
  * Wu et al. (2024) W. Wu, X. Feng, Z. Gao, and Y. Kan.  Smart: Scalable multi-agent real-time simulation via next-token prediction.  _arXiv preprint arXiv:2405.15677_ , 2024.
  * Zhou et al. (2023) Z. Zhou, J. Wang, Y.-H. Li, and Y.-K. Huang.  Query-centric trajectory prediction.  In _CVPR_ , 2023.

## Appendix

| mAP | minADE | minFDE | MissRate
---|---|---|---|---
Vehicle | 0.4745 | 0.7526 | 1.5107 | 0.1489
Pedestrian | 0.4827 | 0.3455 | 0.7251 | 0.0756
Cyclist | 0.3898 | 0.7095 | 1.4299 | 0.1865
Avg | 0.4490 | 0.6025 | 1.2219 | 0.1370

Table 3: Evaluation of our MTR model on the WOMD Motion Prediction dataset
(validation). Figure 5: Three simulated scenarios (top to bottom) at different
timesteps (left to right) showcasing multi-modal behavior of agents. In the
top and bottom simulation, the dark green car takes the left turn. However, in
the middle simulation, it turns right. The green car, in the top and middle
simulation, attempts the lane change to the left as the cars in front wait at
the signal. In the middle simulation, the same green car comes to a stop in
the same lane behind the traffic.  Figure 6: Two simulated scenarios (top to
bottom) at different timesteps (left to right). In the top simulation, the
light green car, attempting to take the right turn, stops and respects the
teal car’s right of way. Whereas, in the bottom simulation, the light green
car, quickly takes the right turn.  Figure 7: Simulated scenario at different
timesteps(left to right). The dark green car, about to take the free right,
waits for the pedestrian to cross the road.  Figure 8: Simulated scenario at
different timesteps(left to right). The dark green car merges to the left lane
as the right lane comes to an end.  Figure 9: Simulated scenario at different
timesteps(left to right). The cars come to stop at the signal. Additionally,
the teal car starts to slow down and stops as the golden car in front is
waiting at the signal.  Figure 10: Simulated scenario at different
timesteps(left to right). The red car waits for the purple car to pass before
taking an unprotected left turn.  Figure 11: Two simulated scenarios (top to
bottom) at different timesteps (left to right). The top simulation shows the
brown car entering the parking lot and driving straight. In the bottom
simulation, the brown car attempts to park besides the pink car. We can also
observe the dark purple car, waiting for the teal car to complete the U-turn
before taking the free right.