#  Shock cooling and breakout emission for optical flares associated with
gravitational wave events

Hiromichi Tagawa11affiliationmark:  , Shigeo S Kimura22affiliationmark:
33affiliationmark:  , Zoltán Haiman44affiliationmark:  55affiliationmark:  ,
Rosalba Perna66affiliationmark:  77affiliationmark:  , Imre
Bartos88affiliationmark:  11affiliationmark:  Shanghai Astronomical
Observatory, Shanghai, 200030, People′s Republic of China
22affiliationmark:  Astronomical Institute, Graduate School of Science, Tohoku
University, Aoba, Sendai 980-8578, Japan
33affiliationmark:  Frontier Research Institute for Interdisciplinary
Sciences, Tohoku University, Sendai 980-8578, Japan
44affiliationmark:  Department of Astronomy, Columbia University, 550 W. 120th
St., New York, NY, 10027, USA
55affiliationmark:  Department of Physics, Columbia University, 550 W. 120th
St., New York, NY, 10027, USA
66affiliationmark:  Department of Physics and Astronomy, Stony Brook
University, Stony Brook, NY 11794-3800, USA
77affiliationmark:  Center for Computational Astrophysics, Flatiron Institute,
New York, NY 10010, USA
88affiliationmark:  Department of Physics, University of Florida, PO Box
118440, Gainesville, FL 32611, USA
E-mail<EMAIL_ADDRESS>
###### Abstract

The astrophysical origin of stellar-mass black hole (BH) mergers discovered
through gravitational waves (GWs) is widely debated. Mergers in the disks of
active galactic nuclei (AGN) represent promising environments for at least a
fraction of these events, with possible observational clues in the GW data. An
additional clue to unveil AGN merger environments is provided by possible
electromagnetic emission from post-merger accreting BHs. Associated with BH
mergers in AGN disks, emission from shocks emerging around jets launched by
accreting merger remnants is expected. In this paper we compute the properties
of the emission produced during breakout and the subsequent adiabatic
expansion phase of the shocks, and we then apply this model to optical flares
suggested to be possibly associated with GW events. We find that the majority
of the reported flares can be explained by the breakout and the shock cooling
emission. If these events are real, then the merging locations of binaries are
constrained depending on the emission processes. If the optical flares are
produced by shock cooling emission, they would display moderate color
evolution, possibly color variations among different events, a positive
correlation between the delay time and the duration of flares, and
accompanying breakout emission in X-ray bands before the optical flares. If
the breakout emission dominates the observed lightcurve, it is expected that
the color is distributed in a narrow range in the optical band, and the delay
time from GW to electromagnetic emission is longer than $\sim 2$ days. Hence,
further explorations of the distributions of delay times, color evolution of
the flares, and associated X-ray emission will be useful to test the proposed
emission model for the observed flares.

###### Subject headings:

transients – stars: black holes –galaxies: active

##  1\. Introduction

The astrophysical pathways to black hole (BH) mergers discovered by the LIGO
(Aasi et al., 2015), Virgo (Acernese et al., 2015), and KAGRA (Akutsu et al.,
2021) gravitational wave (GW) observatories have been actively debated (Barack
et al., 2019). Various scenarios have been proposed, including isolated binary
evolution (e.g. Dominik et al., 2012; Kinugawa et al., 2014; Belczynski et
al., 2016; Spera et al., 2019; Tagawa et al., 2018; Tanikawa et al., 2022)
evolution of triple or quadruple systems (e.g. Silsbee & Tremaine, 2017;
Antonini et al., 2017; Michaely & Perets, 2019; Fragione & Kocsis, 2019;
Martinez et al., 2022; Bartos et al., 2023), dynamical evolution in star
clusters (e.g. Portegies Zwart & McMillan, 2000; Samsing et al., 2014; O’Leary
et al., 2016; Rodriguez et al., 2016; Banerjee, 2017; Kumamoto et al., 2018;
Rasskazov & Kocsis, 2019; Perna et al., 2019; Fragione & Kocsis, 2018;
Antonini et al., 2022), and compact objects in active galactic nucleus (AGN)
disks (e.g. Bartos et al., 2017; Stone et al., 2017; McKernan et al., 2018;
Tagawa et al., 2020b; Arca Sedda et al., 2023).

In an AGN disk, BHs are embedded in due to capture via dynamical interactions
between the nuclear star cluster and the AGN disk (Ostriker, 1983; Miralda-
Escudé & Kollmeier, 2005; Generozov & Perets, 2022; Wang et al., 2023) and by
in-situ star formation (Levin & Beloborodov, 2003; Goodman & Tan, 2004;
Milosavljević & Loeb, 2004; Derdzinski & Mayer, 2023; Chen et al., 2023b). The
AGN disk environment helps to bring the BHs closer together (Bellovary et al.,
2016; Tagawa et al., 2020b; Derdzinski et al., 2019, 2021; Li et al., 2021;
Grishin et al., 2023) and hence form binaries (Li et al., 2022; Boekholt et
al., 2023; Rowan et al., 2022; Li et al., 2023; DeLaurentiis et al., 2023;
Rozner et al., 2023; Whitehead et al., 2023; Rowan et al., 2023; Qian et al.,
2023) which then merge over relatively short time scales. Comparisons to the
observed BH masses, spins and merger rate indicate that a sizable fraction of
the observed mergers may indeed originate in AGN disks (Gayathri et al., 2021;
Tagawa et al., 2021a; Ford & McKernan, 2022). The AGN channel could also
explain some of the peculiar detections, such as those with a high mass
(Tagawa et al., 2021a; Gayathri et al., 2023) and possibly high eccentricity
(Samsing et al. 2020; Tagawa et al. 2021b; Gayathri et al. 2022; Romero-Shaw
et al. 2022 but see Romero-Shaw et al. 2020, 2023).

Due to the gas-rich merger environment, a key signature of the AGN channel is
the possibility of electromagnetic emission accompanying the GW signal from
the merger (Bartos et al., 2017; Stone et al., 2017). To explore this
possibility, electromagnetic follow-up observations have been carried out for
many of the mergers, with nine counterpart candidates suggested so far,
including seven optical flares (Graham et al., 2020, 2023) and two gamma-ray
flares (Connaughton et al., 2016; Bagoly et al., 2016).

Recently, several studies have investigated the electromagnetic emission from
a variety of transients emerging from AGN disks. Many studies (Zhu et al.,
2021c, a; Perna et al., 2021a; Yuan et al., 2021; Wang et al., 2022; Lazzati
et al., 2022; Ray et al., 2023) focused on the radiation from gamma-ray
bursts, while Perna et al. (2021b) and Zhu et al. (2021b) discussed the
electromagnetic signatures expected from accretion-induced collapse of neutron
stars and white dwarfs, respectively. Yang et al. (2021), Xin et al. (2023),
and Prasad et al. (2023) studied the properties of tidal disruption of stars
by BHs, while Grishin et al. (2021) investigated supernova explosions, and
Bartos et al. (2017) and Stone et al. (2017) estimated the electromagnetic
emission produced by thermal radiation and/or outflow from circum-binary disks
in AGN disks.

Several recent studies have also investigated whether transients from BHs
mergingin AGN disks could explain the optical flare, ZTF19abanrhr (Graham et
al., 2020), associated with the BH merger GW190521 detected in GWs. McKernan
et al. (2019) discussed emission from shocks caused by collision between gas
bound to the merger remnant and unbound gas after recoil kicks due to
anisotropic radiation of GWs. Graham et al. (2020) assessed the net luminosity
and timescales for gas accretion induced by recoil kicks. de Mink & King
(2017) considered flares emerging from shocks in a circum-BH disk due to
recoil kicks. Kimura et al. (2021), Wang et al. (2021a, b), and Chen et al.
(2023a), respectively, considered thermal and non-thermal emission from
bubbles and bubble evolution around BHs formed by strong outflows considering
continuous and episodic super-Eddington accretion. Wang et al. (2021a) further
considered emission from shocks emerging due to interactions of Blandford-
Znajek jets (Blandford & Znajek, 1977) launched from accreting BHs to the
broad line regions, Rodríguez-Ramírez et al. (2023) considered free-free and
bound-free emission from gas shocked due to interaction of the jets and AGN
disk gas, and Tagawa et al. (2023a) estimated gamma-ray, neutrino, and cosmic-
ray emission from internal shocks in the jets. Ashton et al. (2020), Palmese
et al. (2021), Calderón Bustillo et al. (2021), and Morton et al. (2023)
estimated the association significance of ZTF19abanrhr to GW190521.

In this paper, we develop an emission model based on the scenario proposed by
Tagawa et al. (2022, 2023b) and discuss whether or not emission based on this
scenario can explain some of the optical flares reported in Graham et al.
(2023). Tagawa et al. (2022) indicated that a Blandford-Znajek jet can be
produced from BHs embedded in AGN disks and investigated its influence on the
AGN disk structure. Due to the high pressure of the shocks emerging around the
jet, a cavity is created around the BHs. Just before the jet breaks out of the
AGN disk, photons in the shocked gas begin to escape. These photons can be
observed as breakout emission (e.g. Nakar & Sari, 2010), whose properties have
been investigated in Tagawa et al. (2023c). The BHs can maintain the jets even
after they break out from the AGN disk, as long as there is leftover circum-BH
disk gas. Once this is depleted, the BHs can no longer power the jets. This is
then followed by an inactive phase which lasts until gas is replenished onto
the BH, and the jet is launched again, with the cycle hence repeating.

In the case that BHs merge in the cavity while accreting (upper panel of Fig.
1), Tagawa et al. (2023b) (hereafter Paper I) predicted that electromagnetic
emission is often produced associated with the BH merger. This is because the
jet direction can be reoriented following a merger ($\S~{}\ref{sec:model}$),
and strong shocks and emission are produced soon after the jet reorientation.
Paper I investigated the properties of breakout emission from a jet head
associated with BH mergers, and found that this model can explain various
properties, including the luminosity, delay time, duration, and color of the
electromagnetic transients ZTF19abanrhr, GW150914-GBM, and LVT151012, as well
as why the transients began brightening only after the merger. Zhu (2023) and
Zhou & Wang (2023) also estimated the neutrino emission from the breakout of
the jets produced associated with BH mergers. In this paper, we additionally
consider the shock cooling emission, which is produced in a subsequent
adiabatic expansion phase of the shocked gas (e.g. Arnett 1980; Morag et al.
2023 for supernovae and Nakar & Piran 2017; Gottlieb et al. 2018; Hamidani &
Ioka 2023 for gamma-ray bursts), in addition to the breakout emission, which
is produced in an early phase of the shocked gas (see Fig. 1 for a schematic
picture). We present the properties of the emission in both phases, and we
then apply this model to the flares reported in Graham et al. (2023). We
further discuss how to test the proposed emission model in future
observations.

The rest of this paper is organized as follows. In § 2, we describe a model
for producing electromagnetic flares associated with GW emission and a way to
constrain physical properties of the flares from the observations using this
model. We present our results in § 3, discuss how to test the model in § 4,
and summarize our conclusions in § 5.

Figure 1.—  A schematic picture of breakout and cooling envelope emission from
shocks produced by the interaction of AGN gas and a jet launched from a BH
accreting gas in an AGN disk. Since the jet is reoriented at the time of the
merger of two BHs ($\S~{}\ref{sec:model}$), electromagnetic emission is also
produced after the merger.

##  2\. Method

First we describe the model itself. We then specialize to discuss how to
derive the model parameters from the observed properties of the flares, that
is the delay time ($t_{\rm delay}$), the duration ($t_{\rm duration}$), the
luminosity of the flare ($L_{\rm obs}$), the merger remnant mass ($m_{\rm
BH}$), the SMBH mass ($M_{\rm SMBH}$), and the AGN luminosity ($L_{\rm AGN}$).
In our analysis, we assume that the shocks produced by collisions between the
jet and the AGN gas are characterized by non-relativistic regimes. This is
because in both the breakout and the shock cooling emission, flares with the
delay time and the duration of $\gtrsim 10~{}{\rm day}$ (Table 1) are usually
characterized by this regime (Paper I, Tables 2 and 3). As possible processes
for explaining the properties of the optical flares, we consider the breakout
emission from the jet head (the breakout emission scenario) and shock cooling
emission from the cocoon (the shock cooling emission scenario). Note that we
use the Shakura-Sunyaev model (Shakura & Sunyaev, 1973) for the accretion disk
in the shock cooling emission scenario, and the Thompson disk model (Thompson
et al., 2005) in the breakout emission scenario. This is because the position
of the BH from the central supermassive BH (SMBH) is sub-parsec for the former
(Table 2) and a few parsec for the latter case (Table 3), and the mechanisms
of angular momentum transfer and the disk properties are likely different in
the two regions, with the Thompson disk model better suited than the Shakura-
Sunyaev one to describe the outer disk.

###  2.1. Shock formation

Shocked gas is responsible for both the breakout and the cooling emission, and
hence we begin by describing the process of shock formation.

We first describe the accretion rate onto BHs in an AGN disk, which is
evaluated by the Bondi-Hoyle-Lyttleton rate. For a BH embedded in a cool AGN
disk, the Bondi-Hoyle-Lyttleton radius ($r_{\rm BHL}$) is large, and usually
exceeds the scale height of the AGN disk ($H_{\rm AGN}$) and the Hill radius
($r_{\rm Hill}$). Accounting for the geometrical limitation of the capture
regions by the shear motion and the vertical height of the AGN disk, the
capture rate of gas by the BH is given by

$\displaystyle{\dot{m}}_{\rm BHL}=$ $\displaystyle f_{\rm c}r_{\rm w}r_{\rm
h}\rho_{\rm AGN}(c_{\rm s,AGN}^{2}+v_{\rm BH}^{2}+v_{\rm sh}^{2})^{1/2}$
$\displaystyle\simeq$ $\displaystyle~{}3\times 10^{-4}\,{\rm M}_{\odot}/{\rm
yr}~{}\left(\frac{f_{\rm c}}{10}\right)\left(\frac{H_{\rm AGN}}{0.003\,{\rm
pc}}\right)\left(\frac{R_{\rm BH}}{1~{}{\rm pc}}\right)^{1/2}$
$\displaystyle\left(\frac{\rho_{\rm AGN}}{4\times 10^{-17}\,{\rm
g/cm^{3}}}\right)\left(\frac{m_{\rm BH}}{10~{}{\rm
M}_{\odot}}\right)^{2/3}\left(\frac{M_{\rm SMBH}}{10^{6}\,{\rm
M}_{\odot}}\right)^{-1/6},$ (1)

(e.g. Tanigawa & Tanaka, 2016; Choksi et al., 2023), where $\rho_{\rm AGN}$ is
the gas density and $c_{\rm s,AGN}$ is the sound speed of the AGN disk at the
position of the BH ($R=R_{\rm BH}$, where $R$ is the distance from the SMBH),
$v_{\rm BH}$ is the velocity of the BH with respect to the local motion of the
AGN disk, $v_{\rm sh}=r_{\rm w}\Omega$ is the shear velocity at the capture
radius $r_{\rm w}={\rm min}(r_{\rm BHL},\,r_{\rm Hill})$, $\Omega=(GM_{\rm
SMBH}/R_{\rm BH}^{3})^{1/2}$ is the angular velocity of the BH, $r_{\rm
h}={\rm min}(r_{\rm w},\,H_{\rm AGN})$ is the capture height, $G$ is the
gravitational constant, and $f_{\rm c}\sim 10$ is a normalization constant
(Tanigawa & Watanabe, 2002). In the second line of Eq. (2.1), we assume
$v_{\rm BH}<c_{\rm s,AGN}$ and $v_{\rm BH}<v_{\rm sh}$. Note that $r_{\rm
w}(c_{\rm s,AGN}^{2}+v_{\rm sh}^{2})^{1/2}\approx r_{\rm Hill}H_{\rm
AGN}\Omega$ is used to derive the right hand side, which is approximately
satisfied regardless of whether $c_{\rm s,AGN}$ is larger or smaller than
$v_{\rm sh}$. By considering the reduction or enhancement with respect to the
Bondi-Hoyle-Lyttleton rate, we parameterize the fraction of the accretion rate
onto the BH (${\dot{m}}$) over the Bondi-Hoyle-Lyttleton rate (${\dot{m}}_{\rm
BHL}$) by $f_{\rm acc}={\dot{m}}/{\dot{m}}_{\rm BHL}$ as in Paper I. For
example, low $f_{\rm acc}$ may be predicted due to winds from an accretion
disk with a super-Eddington rate. On the other hand, recent simulations
suggest that the conversion to wind is moderate (Kitaki et al., 2021) for
accretion flows in which the circularization radius (where gas is circularized
after being captured by a BH) is much larger than the trapping radius (within
which photons are advected to a BH without escaping), as is the case for BHs
embedded in an AGN disk. In addition, the accretion rate onto a BH in a cavity
during the active phases is estimated to be lower by a factor of a few
compared to that without a cavity (Tagawa et al., 2022).

From rapidly accreting and spinning BHs in an AGN disk, a Blandford-Znajek jet
is expected to be launched, as outlined in Appendix A.1 of Tagawa et al.
(2022). The jet kinetic luminosity ($L_{\rm j}$) is proportional to the mass
accretion rate onto the BH (${\dot{m}}$) as

$\displaystyle L_{\rm j}=\eta_{\rm j}{\dot{m}}c^{2},$ (2)

where $c$ is the speed of light, and $\eta_{\rm j}$ is the jet conversion
efficiency, which is approximated by $\eta_{\rm j}\sim a_{\rm BH}^{2}$ for a
magnetically dominated jet (Tchekhovskoy et al., 2010; Narayan et al., 2021),
$a_{\rm BH}$ is the dimensionless spin of the BH, and $a_{\rm BH}\sim 0.7$ for
the merger remnants (The LIGO Scientific Collaboration et al., 2021).

At a BH merger ($t=0$, where $t$ is the time from the merger), the jet
direction is reoriented and can collide with unshocked AGN gas in the
following ways. Once two BHs merge, the BH spin direction is reoriented if the
angular momentum direction of the merging binary is misaligned with respect to
the spin directions of the merging BHs. This is expected for mergers in an AGN
disk due to frequent binary-single interactions (Tagawa et al., 2020a) and/or
inhomogeneity of AGN disks. Since the jet is injected in the direction of the
BH spin, if the jet is not aligned with the circum-BH disk due to a strong jet
power (Polko & McKinney, 2017), the jet propagates in the direction of the BH
spin, and can collide with AGN gas. Even if the jet aligns with the angular
momentum direction of the circum-BH disk due to magnetic interactions on
average, the jet can precess by interacting with magnetic fields while the
angular momentum directions of the BH spin and the circum-BH disk are still
misaligned with one other (Liska et al., 2018, 2021). Due to the precession,
the jet can collide with unshocked gas after merger during the first
precession cycle (after that, the opening angle of the cavity becomes wider
than that of the precession). The other possibility is that once BHs merge, a
merger remnant BH receives a recoil kick in an almost random direction due to
anisotropic radiation of GWs. Then shocks form in a circum-BH disk within
$\lesssim 10^{13}~{}{\rm cm}$ (Tagawa et al., 2023b), and shocked gas accretes
onto the remnant with the angular momentum direction being modified as a
result of shocks (Rossi et al., 2010). Due to magnetic interactions (Liska et
al., 2021), the jet is then aligned with the angular momentum direction of the
circum-BH disk, which is in turn misaligned with respect to the jet direction
before the merger, and the jet can therefore collide with AGN gas.

Once the jet collides with unshocked AGN gas around a BH, two shocks form: a
forward shock propagating in the AGN disk and a reverse shock in the jet. The
region sandwiched by the two shocks is called the jet head. In the jet
collimated regime considered in our work, the dimensionless velocity of
shocked gas in the jet head at the shock breakout is estimated as (Bromberg et
al., 2011)

$\displaystyle\beta_{\rm h}\sim\left(\frac{L_{\rm j}}{\rho_{\rm AGN}t_{\rm
break}^{2}\theta_{0}^{4}c^{5}}\right)^{1/5},$ (3)

where $\theta_{0}$ is the opening angle of the injected jet, and $t_{\rm
break}$ is the delay time between the production of the jet and the shock
breakout, and is roughly given by

$\displaystyle t_{\rm break}$ $\displaystyle\sim\frac{3}{5}\frac{H_{\rm
AGN}}{\beta_{\rm FS}c}f_{\rm corr}$ $\displaystyle\sim 1\,{\rm
yr}\left(\frac{H_{\rm AGN}}{5\times
10^{16}\,\mathrm{cm}}\right)\left(\frac{\beta_{\rm
FS}}{0.1}\right)^{-1}\left(\frac{f_{\rm corr}}{3}\right).$ (4)

Here $\beta_{\rm FS}\sim(7/6)\beta_{\rm h}$ is the dimensionless velocity of
the forward shock of the jet head, and $f_{\rm corr}$ is the correction factor
for the delay time due to the inclination of the jet and the geometry of the
cavity. Ignoring geometrical corrections due to the existence of the cavity,
we assume $f_{\rm corr}={\rm min}[1/{\rm cos}i,1/\theta_{0}{\rm sin}i]$,
considering both the cases in which the shocked gas (cocoon) breaking out of
the AGN disk is from its head and from its sides, where $i$ is the angle
between the jet and the angular momentum direction of the AGN disk. With this
prescription, $f_{\rm corr}$ ranges between $\sim 1$–$1/\theta_{0}$.

After the shock breakout, radiation is produced. Early emission is
characterized by the breakout emission as described in § 2.3 and Paper I,
while later emission is characterized by the shock cooling radiation described
in § 2.2. Since breakout emission is associated with the shock propagation,
both non-thermal and thermal emissions are expected. On the other hand, we
only consider thermal emission from shock cooling, and neglect any non-thermal
emission from expanding ejecta. The latter can also be bright if additional
strong shocks form when the ejecta collide with the interstellar medium. We
use the non-thermal component of the breakout emission and the thermal shock
cooling emission to model the flares reported by Graham et al. (2023), since
these are bright in optical bands. Note that the peak frequency of the thermal
breakout emission from the jet head falls above the X-ray bands, which is
constrained by the duration of the flare (e.g. Eqs. 2.3.1 and 37 below), and
hence it cannot reproduce the optical flares.

###  2.2. Shock cooling emission

####  2.2.1 Physical model

In the shock cooling emission, photons diffuse and are released from deep
inside the shocked gas, as observed in supernovae (Arnett, 1980). To present
the properties of this emission, we first describe the evolution of the
shocked gas.

At the breakout of the shocked gas (cocoon) at $t=t_{\rm break}$, the thermal
energy of the cocoon is roughly given by

$\displaystyle e_{\rm BO}$ $\displaystyle=\frac{1}{2}m_{\rm BO}\beta_{\rm
c}^{2}c^{2},$ (5)

where $\beta_{\rm c}\simeq\beta_{\rm h}\theta_{0}$ is the dimensionless
velocity of the cocoon (Bromberg et al., 2011), and $m_{\rm BO}$ is the mass
of the cocoon at the breakout and is roughly given by

$\displaystyle m_{\rm BO}$ $\displaystyle\simeq 2\pi H_{\rm AGN}^{3}f_{\rm
corr}^{3}\theta_{0}^{2}\rho_{\rm AGN},$ (6)

considering the cylindrical shape of the cocoon with the height of $H_{\rm
AGN}f_{\rm corr}$ and the radius of $H_{\rm AGN}f_{\rm corr}\theta_{0}$. After
the shock passage (even before the shock breakout), the internal energy of the
shocked gas is converted to kinetic energy due to the expansion caused by the
radiation pressure. Afterwards, the shocked ejecta expands nearly spherically
with velocity $v_{\rm ej}\sim\beta_{\rm c}c$. As the ejecta expands to size
$R$, the optical depth of the spherically expanding ejecta declines as (e.g.
Sapir & Waxman, 2017; Nakar & Piran, 2017)

$\displaystyle\tau_{\rm ej}\simeq\frac{\kappa_{\rm ej}m_{\rm BO}}{4\pi
R^{2}},$ (7)

where $\kappa_{\rm ej}$ is the ejecta’s opacity. We adopt the Thomson
scattering opacity of $\kappa_{\rm ej}\sim 0.4\,{\rm g/cm^{2}}$ assuming
ionized gas. If the ejecta is initially optically thick, as assumed in our
fiducial model, photons deep inside the ejecta can be diffused out once the
optical depth is reduced to

$\displaystyle\tau_{\rm ej}=c/v_{\rm ej}$ (8)

at the time

$\displaystyle t_{\rm diff}=\left(\frac{\kappa_{\rm ej}m_{\rm BO}}{4\pi
cv_{\rm ej}}\right)^{1/2},$ (9)

when the corresponding radius is $R_{\rm diff}=t_{\rm diff}v_{\rm ej}$ and the
density

$\displaystyle\rho_{\rm diff}=\frac{c}{\kappa_{\rm ej}R_{\rm diff}v_{\rm
ej}}=\frac{c}{\kappa_{\rm ej}t_{\rm diff}v_{\rm ej}^{2}}.$ (10)

Due to the adiabatic expansion, the thermal energy at the diffusion radius is
reduced by a factor of $\sim(R_{\rm BO}^{3}/R_{\rm diff}^{3})^{\gamma-1}$
where $\gamma=4/3$ is the adiabatic index for the radiation pressure dominated
gas, $R_{\rm BO}=(V_{\rm BO}/4\pi)^{1/3}$ is the typical size, and $V_{\rm
BO}$ is the volume of the cocoon at the shock breakout. Hence, the luminosity
at $\tau_{\rm ej}=c/v_{\rm ej}$ is related to $e_{\rm BO}$ as

$\displaystyle L_{\rm SC}=\frac{e_{\rm BO}}{t_{\rm diff}}\frac{R_{\rm
BO}}{R_{\rm diff}}.$ (11)

Using $R_{\rm BO}$ in Eq. (11), the unperturbed AGN density at the position of
the BH ($R=R_{\rm BH}$) can be estimated via

$\displaystyle\rho_{\rm AGN}\sim\rho_{\rm diff}\left(\frac{R_{\rm
diff}}{R_{\rm BO}}\right)^{3}.$ (12)

Furthermore, from the AGN density, the inflow rate of the AGN disk at
$R=R_{\rm BH}$ can be calculated via

$\displaystyle{\dot{M}}_{\rm inflow}=4\pi H_{\rm AGN}^{3}\rho_{\rm
AGN}\Omega\alpha$ (13)

assuming an alpha viscosity with parameter $\alpha$, and where $\Omega$ is the
orbital angular velocity of the BH around the SMBH. Here, the AGN luminosity
($L_{\rm AGN}$) is related to the inflow rate as

$\displaystyle\frac{L_{\rm AGN}}{\eta_{\rm rad}c^{2}}={\dot{M}}_{\rm
SMBH}=f_{\rm cons}{\dot{M}}_{\rm inflow}$ (14)

where ${\dot{M}}_{\rm SMBH}$ is the accretion rate onto the SMBH, $\eta_{\rm
rad}$ is the radiation efficiency, and $f_{\rm cons}\equiv{\dot{M}}_{\rm
SMBH}/{\dot{M}}_{\rm inflow}\leq 1$ is the fraction of the inflow rate at
$R=R_{\rm BH}$ over the accretion rate onto the SMBH.

The radiation temperature is determined as follows. At $\tau_{\rm ej}=c/v_{\rm
ej}$, the energy density of the radiation within the ejecta shell is

$\displaystyle u_{\gamma}=\frac{L_{\rm SC}\tau_{\rm ej}}{4\pi R_{\rm
diff}^{2}c}.$ (15)

The blackbody temperature of the radiation is then given by

$\displaystyle T_{\rm BB}=(u_{\gamma}/a)^{1/4},$ (16)

where $a$ is the radiation constant. Note that the radiation pressure
dominates the gas pressure at $t\lesssim t_{\rm diff}$ for the models
considered in this paper. Since the flares in Graham et al. (2023) were found
by the Zwicky Transient Facility (ZTF), and the observed optical luminosity
($L_{\rm obs}$) can be roughly estimated as the total observed energy of the
flare in the optical band divided by the duration (provided in Table 3 of
Graham et al. 2023), we pose that $L_{\rm obs}$ needs to match the luminosity
in the ZTF bands as

$\displaystyle L_{\rm obs}=L_{\rm SC}\left.\int^{\nu_{\rm up}}_{\nu_{\rm
down}}\frac{2h\nu^{3}}{c^{2}}\frac{d\nu}{e^{h\nu/k_{\rm B}T_{\rm
BB}}-1}\middle/\sigma T_{\rm BB}^{4},\right.$ (17)

where $\nu_{\rm up}=c/(400~{}{\rm nm})$ and $\nu_{\rm down}=c/(700~{}{\rm
nm})$ are the upper and lower limits for the frequencies observed by the $r$
and $g$ bands of the ZTF. Note that the luminosities in the $r$ and $g$ bands
of the ZTF are not directly derived from Eq. (17). To calculate them (panels c
and d of Fig. 3), we adjust the limits of the integral in this equation to
cover the frequency range of the relevant band.

####  2.2.2 Derivation of model parameters

Using Eqs. (2.1)–(17), we can then determine the 17 quantities $\beta_{\rm
c}$, $\rho_{\rm AGN}$, $R_{\rm BH}$, $H_{\rm AGN}$, $L_{j}$, ${\dot{m}}_{\rm
BHL}$, $t_{\rm break}$, $e_{\rm BO}$, $m_{\rm BO}$, $\tau_{\rm ej}$,
$\rho_{\rm diff}$, $R_{\rm BO}$, $R_{\rm diff}$, $u_{\gamma}$, $T_{\rm BB}$,
$L_{\rm SC}$, and ${\dot{M}}_{\rm inflow}$ given the input parameters,
$\theta_{0}$, $f_{\rm acc}$, $a_{\rm BH}$, $\alpha$, $\eta_{\rm rad}$, $f_{\rm
corr}$, and $f_{\rm cons}$, and the observed properties, $L_{\rm obs}$ $L_{\rm
AGN}$, and $t_{\rm diff}$. For example, Combining Eqs. (2.1)–(14), the
velocity of the cocoon is calculated as

$\displaystyle\beta_{\rm c}=\left[\frac{32\pi^{4}f_{\rm jet/BHL}^{3}f_{\rm
cons}Ac^{10}G^{2}m_{\rm BH}^{2}t_{\rm diff}^{2}\alpha\eta_{\rm rad}}{9L_{\rm
AGN}L_{\rm SC}^{3}\kappa_{\rm ej}^{4}f_{\rm corr}^{6}}\right]^{\delta},$ (18)

where $f_{\rm jet/BHL}\equiv f_{\rm c}f_{\rm acc}\eta_{\rm j}$ is a parameter
related to the jet power, and $A$ and $\delta$ are variables depending on the
velocity of the cocoon. For $\beta_{\rm c}/\theta_{0}<1$,
$A=(35/18)^{6}\theta_{0}^{-3}$ and $\delta=1/2$, and otherwise (as long as the
jet is in the collimated regime) $A=\theta_{0}$ and $\delta=1/8$. To determine
$\beta_{\rm c}$ using Eq. (18), $L_{\rm SC}$ needs to be derived using Eqs.
(15)–(17). To consistently solve these equations, we calculate $\beta_{\rm c}$
in an iterative way using Newton’s method. Then, by incorporating $\beta_{\rm
c}$ into Eqs.(2.1)–(17), the other parameters can also be determined.

We further adjust $f_{\rm jet/BHL}$ so that the scale height of the AGN disk
is equal to the height expected for the Shakura-Sunyaev disk, given $\alpha$.
Such disk structure is expected to be realized in regimes where the Toomre
parameter (Eq. 20) becomes greater than $1$.

###  2.3. Breakout emission

####  2.3.1 Physical model

Long before photons deep inside the ejecta escape ($t\sim t_{\rm break}+t_{\rm
diff}$), at around the time that the shock arrives to the surface of the
ejecta ($t\sim t_{\rm break}$), bright breakout emission is released. Since
breakout emission of the jet head is associated with the shock propagation,
both non-thermal and thermal emission are expected. Such non-thermal emission
can explain the optical flares found by Graham et al. (2020), as investigated
in Paper I. In the following we describe how to reconstruct the model
properties for the breakout emission.

Photons inside the shock start to diffuse out from the AGN disk and the
breakout emission begins to be released when the photon diffusion time from
the shock front, $t_{\rm diff,sh}\sim~{}d_{\rm edge}^{2}\kappa_{\rm
sh}\rho_{\rm AGN}/c$, becomes equal to the dynamical timescale of the shock,
$t_{\rm dyn,edge}\sim~{}d_{\rm edge}/\beta_{\rm FS}c$, where $d_{\rm edge}$ is
the thickness of the AGN disk above the shock, and $\kappa_{\rm sh}$ is
opacity of the shocked gas. By equating these timescales, the thickness at the
breakout is given by $d_{\rm edge,BO}\sim~{}1/(\beta_{\rm FS}\kappa_{\rm
sh}\rho_{\rm AGN})$, and the duration of the emission from a breakout shell is

$\displaystyle t_{\rm diff,BO}$ $\displaystyle\sim 1/(\beta_{\rm
FS}^{2}c\kappa_{\rm sh}\rho_{\rm AGN})$ $\displaystyle\sim 3\,{\rm
yr}\left(\frac{\beta_{\rm FS}}{0.1}\right)^{-2}\left(\frac{\rho_{\rm
AGN}}{1\times 10^{-16}\,{\rm g~{}cm^{-3}}}\right)^{-1},$ (19)

where we adopt $\kappa_{\rm sh}\sim 0.4\,{\rm g/cm^{2}}$ considering the
ionization of gas by photons released from the shocks (Nakar & Sari, 2010).

If we assume that the AGN disk is gravitationally unstable at the position of
the BH, which is expected for the delay time and the duration of the flares
(Paper I), the density of the AGN disk is related to the position of the BH
through (e.g. Thompson et al., 2005)

$\displaystyle\frac{\Omega^{2}}{\sqrt{2}\pi G\rho_{\rm AGN}}=Q\simeq 1,$ (20)

where $Q$ is the Toomre parameter. Note that AGN disks at the BH positions are
predicted to be Toomre unstable (Table 2, e.g. Haiman et al., 2009).

####  2.3.2 Derivation of model parameters

Using Eqs.(2.1)–(2.1) and (2.3.1)–(20), we can determine the 6 variables
$\beta_{\rm h}$, $\rho_{\rm AGN}$, $R_{\rm BH}$, $H_{\rm AGN}$, $L_{j}$, and
${\dot{m}}_{\rm BHL}$, given $\theta_{0}$, $f_{\rm acc}$, and $\eta_{\rm j}$.
For example, the velocity of the jet head is calculated by

$\displaystyle\beta_{\rm h}=\left[\frac{5(7/6)^{4/3}f_{\rm jet/BHL}m_{\rm
BH}^{2/3}G^{1/2}\left(\frac{c\kappa_{\rm sh}t_{\rm
diff,BO}}{\sqrt{2}\pi}\right)^{1/6}}{3^{5/3}f_{\rm corr}t_{\rm
break}\theta_{0}^{4}c^{2}}\right]^{3/11}.$ (21)

Then, by incorporating $\beta_{\rm h}$ into Eqs.(2.1)–(2.1) and (2.3.1)–(20),
the other parameters can be also determined.

In the case of breakout emission, $R_{\rm BH}$ is found to be large (§ 3.2).
At such large scales, efficient transfer of angular momentum of the AGN gas is
required for SMBH accretion (Thompson et al., 2005). Following Thompson et al.
(2005), we assume that the inflow rate of the AGN disk is parameterized by

$\displaystyle{\dot{M}}_{\rm inflow}=4\pi R_{\rm BH}H_{\rm AGN}^{2}\rho_{\rm
AGN}\Omega m.$ (22)

###  2.4. Evolution

In this section, we describe the evolution of the luminosity and the
temperature for thermal emission from the cocoon in the breakout and shock
cooling emission phases. Referring to previous studies (Arnett, 1980; Nakar &
Sari, 2010; Piro et al., 2021; Morag et al., 2023), we assume that the
luminosity evolves as

$\displaystyle L(t^{\prime})\sim\frac{e_{\rm BO}}{t_{\rm dyn}}$
$\displaystyle\left\\{\begin{array}[]{l}0~{}~{}~{}~{}\mathrm{for}~{}~{}t^{\prime}\lesssim
0,\\\ 1~{}~{}~{}~{}\mathrm{for}~{}~{}0\lesssim t^{\prime}\leq t_{\rm pl},\\\
\left(\frac{t^{\prime}}{t_{\rm
pl}}\right)^{-4/3}~{}~{}~{}~{}\mathrm{for}~{}~{}t_{\rm pl}\leq
t^{\prime}\lesssim t_{\rm sph},\\\ \left(\frac{t_{\rm sph}}{t_{\rm
pl}}\right)^{-4/3}\left(\frac{t^{\prime}}{t_{\rm
sph}}\right)^{-\frac{2.28n-2}{3(1.19n+1)}}\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathrm{for}~{}~{}t_{\rm sph}\leq
t^{\prime}\leq t_{\rm diff},\\\ \frac{t_{\rm dyn}}{t_{\rm diff}}\frac{R_{\rm
BO}}{R_{\rm diff}}{\rm exp}\left[-\frac{1}{2}\left(\frac{t^{\prime 2}}{t_{\rm
diff}^{2}}-1\right)\right]\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathrm{for}~{}~{}t_{\rm diff}\leq
t^{\prime},\end{array}\right.$ (30)

where $t^{\prime}=t-t_{\rm break}$, $t_{\rm pl}=1/(\beta_{\rm
c}^{2}c\kappa_{\rm ej}\rho_{\rm AGN})$ is the duration of emission from a
breakout shell in the cocoon, $t_{\rm sph}\simeq R_{\rm BO}/\beta_{\rm c}c$ is
the transition between the planar and spherical geometries of the breakout
shell, $t_{\rm dyn}=H_{\rm AGN}/\beta_{\rm c}c$ is the dynamical timescale,
and $n$ is the power law slope of the vertical AGN gas density profile at the
height at which photons begin to break out. The second and third rows
correspond to the luminosity in the planar phase (before the shocked gas
doubles its radius by expansion) for the breakout emission, the fourth row
corresponds to that in a spherical phase (after the shocked gas doubles and
before the shock cooling emission phase), and the fifth row corresponds to
that in a shock cooling emission phase.

Similarly following previous studies, we assume that the temperature evolves
as

$\displaystyle T(t^{\prime})\sim T_{\rm BB,0}$
$\displaystyle\left\\{\begin{array}[]{l}1~{}~{}~{}~{}\mathrm{for}~{}~{}0\lesssim
t^{\prime}\leq t_{\rm pl},\\\ \left(\frac{t^{\prime}}{t_{\rm
pl}}\right)^{-\frac{2}{3}\frac{9n+5}{17n+9}}~{}~{}~{}~{}\mathrm{for}~{}~{}t_{\rm
pl}\leq t^{\prime}\lesssim t_{\rm sph},\\\ \left(\frac{t_{\rm sph}}{t_{\rm
pl}}\right)^{-\frac{2}{3}\frac{9n+5}{17n+9}}\left(\frac{t^{\prime}}{t_{\rm
sph}}\right)^{n_{\rm T}}\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathrm{for}~{}~{}t_{\rm sph}\leq
t^{\prime}\leq t_{\rm diff},\\\ \frac{T_{\rm BB}}{T_{\rm
BB,0}}\left(\frac{t^{\prime 2}}{t_{\rm
diff}^{2}}-1\right)^{-1/2}~{}~{}~{}~{}\mathrm{for}~{}~{}t_{\rm diff}\leq
t^{\prime},\end{array}\right.$ (36)

where

$\displaystyle T_{\rm BB,0}=(7\rho_{\rm AGN}\beta_{\rm c}^{2}c^{2}/2a)^{1/4}$
(37)

is the black body temperature of the breakout shell, $n_{\rm T}$ is the power-
law index for the temperature evolution in the spherical phase, and $n_{\rm
T}=-(18.48n^{2}+20.69n+6)/[(1.19n+1)(22.32n+17)]$ for an expanding spherical
ejecta (Nakar & Sari, 2010). The validity of these analytical formulae has
been tested by recent numerical simulations (e.g. Piro et al., 2021; Morag et
al., 2023).

Table 1The properties of the ZTF flares suggested to be associated with the GW events. In each pair, the flare is located within the $90\%$ credible volume of and detected within the delay time of $200$ days from the GW event. Here, $M_{\rm SMBH}$ is the SMBH mass, $m_{\rm BH}$ is the merged remnant mass, $t_{\rm duration}$ is the duration of the flare, $t_{\rm delay}$ is the delay time between the detection of the GWs to and the flare, $L_{\rm op}$ is the optical luminosity of the flare, and $L_{\rm AGN}$ is the Eddington ratio of the host AGN. For $\rm J053408.41+085450.6$, the SMBH mass is not constrained and a fiducial mass of $10^{8}~{}{\rm M}_{\odot}$ is used as in Graham et al. (2023).  Pair | LIGO/Virgo alart ID | Name of AGN | ${\rm log_{10}}\left(\frac{M_{\rm SMBH}}{{\rm M}_{\odot}}\right)$ | $m_{\rm BH}~{}[{\rm M}_{\odot}]$ | $t_{\rm duration}~{}[\rm day]$ | $t_{\rm delay}~{}[{\rm day}]$ | $L_{\rm op}~{}[\rm erg/s]$ | $L_{\rm AGN}/L_{\rm Edd}$
---|---|---|---|---|---|---|---|---
1 | $\rm GW190403\\_051519$ | $\rm J124942.30+344928.9$ | 8.6 | 111 | 45.3 | 89 | 44.9 | 0.06
2 | $\rm GW190514\\_065416$ | J124942.30+344928.9 | 8.6 | 68 | 45.3 | 48 | 44.9 | 0.06
3 | $\rm GW190521$ | J124942.30+344928.9 | 8.6 | 164 | 45.3 | 41 | 44.9 | 0.06
4 | $\rm GW190403\\_051519$ | J183412.42+365655.3 | 9.1 | 111 | 41.0 | 110 | 43.8 | 0.02
5 | $\rm GW190424\\_180648$ | J181719.94+541910.0 | 8.0 | 73 | 35.6 | 189 | 44.8 | 0.05
6 | $\rm GW190514\\_065416$ | J224333.95+760619.2 | 8.8 | 68 | 18.3 | 154 | 45.2 | 0.06
7 | $\rm GW190731\\_140936$ | J053408.41+085450.6 | (8.0) | 70.1 | 27.4 | 186 | 44.5 | 0.04
8 | $\rm GW190803\\_022701$ | J053408.41+085450.6 | (8.0) | 65 | 27.4 | 183 | 44.5 | 0.04
9 | $\rm GW190803\\_022701$ | J120437.98+500024.0 | 8.0 | 65 | 47.4 | 191 | 44.7 | 0.02
10 | $\rm GW190909\\_114149$ | J120437.98+500024.0 | 8.0 | 75 | 47.4 | 154 | 44.7 | 0.02
11 | $\rm GW200216\\_220804$ | J154342.46+461233.4 | 9.3 | 81 | 123.4 | 73 | 44.3 | 0.10
12 | $\rm GW200220\\_124850$ | J154342.46+461233.4 | 9.3 | 67 | 123.4 | 69 | 44.3 | 0.10

Since Eqs. (30) and (36) are formulae for an expanding ejecta with a
spherically symmetric evolution, some modifications are required in the cases
of emission from an expanding ejecta with a cylindrical shape like a cocoon,
especially at around the transition between the planar and spherical phases.
For simplicity, we determine $n$ and $n_{\rm T}$ so that the luminosity and
the temperature smoothly evolve at $t^{\prime}=t_{\rm diff}$, respectively.

###  2.5. Parameters

In this section we describe the fiducial values for the model parameters, and
the observed properties of the flares used in the modellings.

As fiducial values, we set the opening angle of the injected jet to
$\theta_{0}=0.2$ (e.g. Hada et al., 2013, 2018; Berger, 2014), the radiation
efficiency to $\eta_{\rm rad}=0.1$, the correction factor for the delay time
to $f_{\rm corr}=3$ (which roughly corresponds to the median value considering
the cavity with the aspect ratio of $\sim 1$ and isotropic jet directions),
the consumption fraction of the inflow rate to $f_{\rm cons}=1$, the alpha-
viscous parameter to $\alpha=0.1$, and the angular momentum transfer parameter
in the outer regions of the AGN disk to $m=0.1$. For the computation of the
shock cooling emission we adjust the value of $f_{\rm jet/BHL}\equiv f_{\rm
c}f_{\rm acc}\eta_{\rm j}$ so that the scale height of the AGN disk is equal
to the height expected for the Shakura-Sunyaev disk (§ 2.2), while for the
breakout emission we set it to $f_{\rm jet/BHL}=10$ assuming $f_{\rm c}=10$
(Tanigawa & Watanabe, 2002), $\eta_{\rm j}=0.5$, and $f_{\rm acc}=2$
considering moderate enhancement of accretion due to shocks caused by recoil
kicks (Appendix B, Tagawa et al., 2023b).

When modeling the flares with breakout emission, we assume that the flare
duration ($t_{\rm duration}$) corresponds to the exponential decay time
($t_{e}$) in Graham et al. (2023). To derive the delay time ($t_{\rm delay}$),
we identify the day at which the flare luminosity peaks using
digitizer111https://automeris.io/WebPlotDigitizer/. Due to the difficulty of
identifying the peak, $t_{\rm delay}$ contains uncertainties. We adopt $M_{\rm
SMBH}$ and $m_{\rm BH}$ from Tables 3 and 4 of Graham et al. (2023). The
optical luminosity is derived by means of the total observed energy in the
optical band divided by the sum of the rise time ($t_{g}$) and the decay time
($t_{e}$) presented in Graham et al. (2023). We calculate the luminosity of
the AGNs ($L_{\rm AGN}$) by inferring the flux of the AGNs at $\sim 4000$Å
from Fig. 6 of Graham et al. (2023), and estimating the luminosity distance
from the redshift of the GW events assuming a value for the Hubble constant of
$67.8~{}{\rm km/s}~{}{\rm s}^{-1}~{}{\rm Mpc}^{-1}$, for the matter density
today of 0.24, and for the cosmological constant today of 0.74 (Planck
Collaboration et al., 2016), and adopting the bolometric correction factor of
$5$ (Duras et al., 2020). The values for the observed quantities adopted in
this paper are listed in Table 1. Note that the delay time ($t_{\rm delay}$)
corresponds to $t_{\rm break}$, and the duration of a flare ($t_{\rm
duration}$) corresponds to $t_{\rm diff,BO}$.

Conversely, when modeling the flares with shock cooling emission, we assume
that the duration of a flare ($t_{\rm duration}$) corresponds to $t_{\rm
diff}$. The delay time between a GW event and an optical flare ($t_{\rm
delay}$) is also on the order of $\sim t_{\rm diff}$ given $t_{\rm break}\ll
t_{\rm diff}$, while we do not use the delay time ($t_{\rm delay}$) to
constrain the model parameters. Note that $t_{\rm delay}\sim t_{\rm duration}$
expected in this scenario is roughly consistent with the properties of the
observed flares (table 1). We assume that the observed optical luminosity of
the flare ($L_{\rm op}$) corresponds to the luminosity in the ZTF bands
($L_{\rm obs}$).

We note that in our model, the physical properties are uniquely determined.
This is because we fixed several input parameters as detailed above, and we do
not directly use the observational data points for parameter fitting. If
instead we allowed the input parameters to vary, we would introduce
degeneracies among several input parameters. The variations of the values of
the fixed input parameters can affect physical properties of the model
especially in the breakout-emission scenario, while less affect them in the
shock cooling emission scenario (see § 4.1 and § 4.2). We believe that such
simple prescriptions are useful to understand typical properties expected from
the model and consider possible tests of the model below.

Table 2 The inferred model parameters for the shock cooling emission scenario. $R_{\rm BH}$ is the distance from the SMBH to the BH, $\beta_{\rm c}$ is the expansion velocity of the cocoon, $h_{\rm AGN}\equiv H_{\rm AGN}/R_{\rm BH}$ and $\rho_{\rm AGN}$ are the aspect ratio and the density of the AGN disk at the position of the BH, respectively, $T_{\rm BB}$ and $\lambda_{\rm peak}$ are the radiation temperature and wavelength of the shock cooling emission, respectively, $L_{\rm jet}$ is the jet power, $f_{\rm jet/BHL}$ is a parameter related to the jet power, $t_{\rm break}$ is the breakout timescale of the jet, and $t_{\rm diff,CBO}$ is the diffusion timescale for the breakout emission from the cocoon.  Pair | $R_{\rm BH}~{}[{\rm pc}]$ | $R_{\rm BH}~{}[R_{\rm g}]$ | $\beta_{\rm c}$ | $h_{\rm AGN}$ | $\rho_{\rm AGN}~{}[\rm g/cm^{3}]$ | $T_{\rm BB}~{}[{\rm K}]$ | $\lambda_{\rm peak}~{}[{\rm nm}]$ | $L_{\rm jet}~{}[{\rm erg/s}]$ | $f_{\rm jet/BHL}$ | $t_{\rm break}~{}[\rm day]$ |  $t_{\rm diff,CBO}~{}[\rm s]$ (Eq. 38)
---|---|---|---|---|---|---|---|---|---|---|---
1 | 0.1 | $5\times 10^{3}$ | 0.2 | $0.004$ | 1$\times 10^{-11}$ | $1.4\times 10^{4}$ | 350 | $1\times 10^{48}$ | 40 | 0.9 | 200
2 | 0.1 | $5\times 10^{3}$ | 0.2 | 0.004 | 1$\times 10^{-11}$ | $1.4\times 10^{4}$ | 350 | $1\times 10^{48}$ | 50 | 0.9 | 200
3 | 0.1 | $5\times 10^{3}$ | 0.2 | 0.004 | 1$\times 10^{-11}$ | $1.4\times 10^{4}$ | 350 | $1\times 10^{48}$ | 30 | 0.9 | 300
4 | 0.08 | $1\times 10^{3}$ | 0.08 | 0.002 | 7$\times 10^{-11}$ | $1.4\times 10^{4}$ | 350 | $1\times 10^{47}$ | 2 | 0.8 | 200
5 | 0.1 | $3\times 10^{4}$ | 0.2 | 0.005 | 1$\times 10^{-12}$ | $1.8\times 10^{4}$ | 270 | $2\times 10^{47}$ | 40 | 2 | $3\times 10^{3}$
6 | 0.06 | $2\times 10^{3}$ | 0.4 | 0.003 | 9$\times 10^{-11}$ | $1.3\times 10^{4}$ | 380 | $2\times 10^{50}$ | 3000 | 0.2 | 7
7 | 0.1 | $2\times 10^{4}$ | 0.1 | 0.005 | 2$\times 10^{-12}$ | $1.9\times 10^{4}$ | 250 | $2\times 10^{47}$ | 30 | 1 | $2\times 10^{3}$
8 | 0.1 | $2\times 10^{4}$ | 0.1 | 0.005 | 2$\times 10^{-12}$ | $1.9\times 10^{4}$ | 250 | $2\times 10^{47}$ | 30 | 1 | $2\times 10^{3}$
9 | 0.3 | $6\times 10^{4}$ | 0.1 | 0.005 | 2$\times 10^{-13}$ | $2.0\times 10^{4}$ | 240 | $5\times 10^{46}$ | 20 | 6 | $5\times 10^{4}$
10 | 0.3 | $6\times 10^{4}$ | 0.1 | 0.005 | 2$\times 10^{-13}$ | $2.0\times 10^{4}$ | 240 | $5\times 10^{46}$ | 20 | 6 | $5\times 10^{4}$
11 | 0.09 | $9\times 10^{2}$ | 0.09 | 0.003 | 2$\times 10^{-10}$ | $9.2\times 10^{3}$ | 520 | $8\times 10^{47}$ | 5 | 1 | 60
12 | 0.09 | $9\times 10^{2}$ | 0.09 | 0.003 | 2$\times 10^{-10}$ | $9.2\times 10^{3}$ | 520 | $8\times 10^{47}$ | 5 | 1 | 50
Table 3 The inferred model parameters for the breakout emission scenario. ${\dot{M}}_{\rm inflow}/{\dot{M}}_{\rm Edd}$ is the Eddington ratio for the gas inflow at the position of the BH, $f_{\rm op/kin}$ is the fraction of the optical luminosity of the flare over the kinetic power of the jet, $f_{\rm BH/NSC}$ is the fraction of $R_{\rm BH}$ over the size of the nuclear star cluster, and $h_{\rm AGN}/h_{\rm AGN,TQM}$ is the ratio of the aspect ratio of the AGN disk derived from our model to that derived assuming the model in Thompson et al. (2005). In the last column, the value of the correction factor for the delay time, at which $h_{\rm AGN}=h_{\rm AGN,TQM}$ is satisfied, is listed.  Pair | $R_{\rm BH}~{}[{\rm pc}]$ | $\beta_{\rm h}$ | $h_{\rm AGN}$ | $\rho_{\rm AGN}~{}[\rm g/cm^{3}]$ | ${\dot{M}}_{\rm inflow}/{\dot{M}}_{\rm Edd}$ | $f_{\rm op/kin}$ | $f_{\rm BH/NSC}$ | $h_{\rm AGN}/h_{\rm AGN,TQM}$ |  $f_{\rm corr}$ for $h_{\rm AGN}=h_{\rm AGN,TQM}$
---|---|---|---|---|---|---|---|---|---
1 | 3.3 | 0.31 | 0.004 | $2\times 10^{-16}$ | 0.06 | 0.36 | 0.38 | 1.05 | 3.1
2 | 3.5 | 0.33 | 0.002 | $1\times 10^{-16}$ | 0.02 | 0.98 | 0.41 | 1.6 | 4.3
3 | 4.0 | 0.41 | 0.002 | $1\times 10^{-16}$ | 0.01 | 0.73 | 0.46 | 1.5 | 4.2
4 | 4.5 | 0.29 | 0.004 | $2\times 10^{-16}$ | 0.05 | 0.02 | 0.44 | 1.03 | 3.1
5 | 1.6 | 0.23 | 0.015 | $4\times 10^{-16}$ | 0.9 | 0.11 | 0.23 | 1.4 | 4.0
6 | 2.4 | 0.23 | 0.008 | $7\times 10^{-16}$ | 0.4 | 0.22 | 0.26 | 1.6 | 4.4
7 | 1.4 | 0.23 | 0.016 | $5\times 10^{-16}$ | 1.2 | 0.05 | 0.21 | 1.6 | 4.3
8 | 1.4 | 0.22 | 0.016 | $5\times 10^{-16}$ | 1.2 | 0.05 | 0.21 | 1.6 | 4.3
9 | 1.7 | 0.23 | 0.014 | $3\times 10^{-16}$ | 0.7 | 0.14 | 0.25 | 1.3 | 3.7
10 | 1.8 | 0.25 | 0.011 | $3\times 10^{-16}$ | 0.4 | 0.16 | 0.26 | 1.1 | 3.3
11 | 8.1 | 0.32 | 0.002 | $6\times 10^{-17}$ | 0.004 | 0.34 | 0.74 | 2.0 | 5.2
12 | 8.0 | 0.31 | 0.001 | $6\times 10^{-17}$ | 0.004 | 0.41 | 0.73 | 1.9 | 4.9
Figure 2.—  The SEDs for shock cooling (orange) and preceding breakout
emission (blue) from a cocoon in the shock cooling emission scenario for pairs
1–12 (left panel), and those for thermal breakout emission (brown) from a jet
head in the breakout emission scenario (right panel) assuming the luminosity
distance of $d_{\rm L}=3~{}{\rm Gpc}$. The spectral energy distributions
(SEDs) for the non-thermal emission is drawn for the model adopted in Fig. 5
of Paper I as a representative example. The dotted blue, purple, red, green,
dark green, and gray lines mark the sensitivity of the ZTF, HiZ-GUNDAM and
Einstein Probe, Chandra and XMM-Newton, NuSTAR, ULTRASAT, and Swift BAT,
respectively, as adopted in Paper I. The vertical cyan lines present the $g$
and $r$ bands of ZTF. Note that the shock cooling emission is too dim ($\ll
10^{-17}~{}{\rm erg/s/cm^{2}}$) in the breakout emission scenario to be
visible in the right panel.  Figure 3.—  The evolution of the luminosity
(panel a), temperature (panel b), apparent magnitude in the $g$ ($m_{g}$,
blue) and $r$ ($m_{r}$, red) bands, AB magnitude at a near-ultraviolet band
(260$~{}{\rm nm}$, black, panel c), and the color ($m_{g}-m_{r}$, panel d) for
thermal emission in the shock cooling emission model for pairs 1 (solid), 5
(dashed), and 11 (dotted). We present pairs 1, 5, and 11 as a fiducial, a
smaller $m_{\rm BH}$, and a longer $t_{\rm duration}$ case, respectively,
among which evolution timescales are derived to be different. Lines are drawn
until $\tau_{\rm ej}=1$ after which our model is invalid.

##  3\. Results

###  3.1. Shock cooling emission

Table 1 shows the list of possible pairs of associations between the GW events
and the electromagnetic flares reported in Graham et al. (2023), along with
their observed properties. Seven electromagnetic flares and twelve pairs for
the possible associations are reported. The number of the pairs (12) is larger
than that of the flares (7), because some flares can be associated with more
than one GW events. Hence, at least five pairs are false associations. Note
that in this paper we do not analyze the two gamma-ray flares possibly
associated with GW events due to their significantly different properties
(Connaughton et al., 2016; Bagoly et al., 2016), while they can be explained
by the thermal breakout emission from the jet head in relativistic regimes
(Tagawa et al., 2023b).

Table 2 shows the distribution of the model parameters when the properties of
the observed flares (Table 1) are modeled with shock cooling emission.

The parameter $f_{\rm jet/BHL}$ is widely distributed depending on the pairs.
$f_{\rm jet/BHL}\sim 1$–$10$ in pairs 4, 11, and 12 is roughly expected for
the Bondi-Hoyle-Lyttleton accretion, while $f_{\rm jet/BHL}\sim 10$–$60$ in
pairs 1–3, 5, and 7–10 can be realized by the enhancement of accretion by
shocks due to recoil kicks at merger (Paper I). However, $f_{\rm jet/BHL}\sim
3000$ in pair 6 is difficult to be realized, since only up to $f_{\rm
jet/BHL}\lesssim 300$ is found to be feasible by considering the radial
surface density profile of a circum-BH disk and appealing to the results of
numerical simulations (see Appendix B of Paper I). On the other hand, even for
pair 6, if we adopt $f_{\rm corr}\sim 1$, $f_{\rm jet/BHL}$ is reduced to
$\sim 300$. This is presumably because both a high accretion rate and a low
inclination have similar influences on the breakout velocity, and most
parameters are mainly characterized by the magnitude of the breakout velocity.
Due to the variation of $f_{\rm corr}$ (between $\sim 1$–$1/\theta_{0}$) and
the degeneracy between $f_{\rm corr}$ and $f_{\rm jet/BHL}$, the cooling model
is poorly constrained or tested by the value of $f_{\rm jet/BHL}$.

The parameter $R_{\rm BH}$ ranges between $0.06$–$0.3$ pc, corresponding to
$900$–$6\times 10^{4}~{}R_{g}$, where $R_{\rm g}=GM_{\rm SMBH}/c^{2}$ is the
gravitational radius of the SMBH. These locations roughly correspond to
migration traps, gap forming regions, and/or slow migration regions, where BH
mergers are predicted to be frequent (Bellovary et al., 2016; Tagawa et al.,
2020b).

The orange lines in the left panel of Fig. 2 show the SEDs for pairs 1–12. The
radiation temperature ranges between $T_{\rm BB}=9200$–$20000~{}{\rm K}$,
corresponding to the peak wavelength $\lambda_{\rm peak}=ch/(3k_{\rm B}T_{\rm
BB})=240$–$520~{}{\rm nm}$, where $h$ is the Planck constant and $k_{\rm B}$
is the Boltzmann constant. Also, the various colors observed in the optical
flares (Fig. 3 of Graham et al. 2023) can be reproduced depending on the
variation of the radiation temperature with changing some input parameters,
such as the direction of the jet ($i$). For example, if we set $i=0^{\circ}$
(the jet being perpendicular to the AGN-disk plane) by using derived values
for $\rho_{\rm AGN}$, $R_{\rm BH}$, and $H_{\rm AGN}$ and the observed
parameters ($m_{\rm BH}$ and $M_{\rm SMBH}$), the radiation temperature of
flares is enhanced by a factor of $\sim 1.7$. This is because for
$i=0^{\circ}$ the shocked mass becomes low, and then photons can escape from
an earlier phase when the radiation temperature is higher. Thus, the variety
of colors can be reproduced if there are temperature variations in the thermal
emission. For synchrotron emission (whose possible contribution to the flares
is discussed in § 3.2), the color is related to the power-law slope of the
injected electrons accelerated by the first-order Fermi process ($p$) as $\nu
L_{\nu}\propto\nu^{(-p+2)/2}$, and $p$ is presumably distributed in a narrow
range for the same phenomena (Panaitescu & Kumar, 2001). On the other hand,
various values for the color can be realized if the breaks of the power laws
in the SEDs coincidentally fall at around the ZTF bands (see the black line in
the right panel of Fig. 2). Note that the distribution of the colors is
currently uncertain due to the shift in the baseline of the flux and the
contribution from the background AGN emission in Fig. 3 of Graham et al.
(2023). If the color is actually distributed in a wide range, the flares are
then easier to model by thermal shock cooling emission rather than by non-
thermal breakout emission. In addition, the model can be also discriminated by
constraining the spectral shape of the emission. This can be done by
simultaneously observing the flares by the ZTF and the Ultraviolet Transient
Astronomy Satellite (ULTRASAT, Shvartzvald et al. 2023, Fig. 2) in the future.

We also predict breakout emission from the cocoon preceding the shock cooling
emission (this is the emission from the cocoon, which is different from the
emission from the jet head discussed in § 3.2). The duration of the breakout
emission from the cocoon,

$\displaystyle t_{\rm diff,CBO}=\frac{1}{\beta_{\rm c}^{2}c\kappa_{\rm
ej}\rho_{\rm AGN}},$ (38)

is distributed in the range $t_{\rm diff,CBO}\sim 10$–$10^{5}~{}{\rm s}$ (the
right most column of Table 2), and the temperature (Eq. 37) is in a range
$2\times 10^{5}$–$2\times 10^{6}~{}{\rm K}$ (blue lines in the right panel of
Fig. 2). Although the duration of the breakout emission is shorter than than
that of the shock cooling emission, it can be detected by future X-ray
surveys, such as HiZ-GUNDAM (Yonetoku et al., 2020) and/or the Einstein Probe
(Yuan & Narayan, 2014).

Fig. 3 a and b show the evolution of the luminosity and the temperature of
emission for pairs 1, 5, and 11. We chose pair 1 as a fiducial example, pair 5
as an event with small $m_{\rm BH}$, and pair 11 as the one with long $t_{\rm
duration}$. The lines are drawn until the phase at which $\tau_{\rm ej}=1$ is
satisfied, since our model for thermal emission is not valid when $\tau_{\rm
ej}<1$.

For the smaller $m_{\rm BH}$ case (dashed lines in Fig. 3), the timescale for
the breakout emission of the cocoon ($t_{\rm pl}$) is longer. This is because
$L_{\rm j}$ is lower due to the lower Bondi-Hoyle-Lyttleton rate (Eqs. 2.1,
2), and $\beta_{\rm c}$ tends to be lower for lower $L_{\rm j}$ (Eqs. 3, 18).
To reproduce $t_{\rm diff}$, $\rho_{\rm AGN}$ needs to be lower for lower
$\beta_{\rm c}$ (Eq. 9). Due to the low $\beta_{\rm c}$ and $\rho_{\rm AGN}$
(Table 2), the breakout timescale of the cocoon ($t_{\rm pl}\propto\beta_{\rm
c}^{-2}\rho_{\rm AGN}^{-1}$) becomes longer. In such case, since photons
diffuse out from the shocked gas earlier, adiabatic expansion is inefficient,
leading to higher temperature.

For the longer $t_{\rm duration}$ model (dotted lines in Fig. 3), the
evolution of the luminosity is slower in the shock cooling regime
($t^{\prime}\gtrsim t_{\rm diff}$).

###  3.2. Breakout emission

Table 3 shows the parameters of the breakout emission model reproducing the
properties of the observed flares. In this model, the parameters are
distributed in narrow ranges. The mergers are predicted to occur at $R_{\rm
BH}\sim 1$–$10~{}{\rm pc}$, the shock velocity ranges within $\beta_{\rm
h}\sim 0.3$–$0.4$, the disk aspect ratio at $R_{\rm BH}$ ranges within $h_{\rm
AGN}\equiv H_{\rm AGN}/R_{\rm BH}\sim 0.001$–$0.02$, the fraction of the
optical luminosity to the kinematic power of the shock is around $L_{\rm
op}/L_{\rm j}\sim 0.02$–$1$, the fraction of the distance from the SMBH to the
size of the nuclear star cluster ($R_{\rm NSC}$) ranges within $R_{\rm
BH}/R_{\rm NSC}\sim 0.2-0.7$ assuming the empirical relation of $R_{\rm
NSC}=8.3~{}{\rm pc}(M_{\rm SMBH}/3.1\times 10^{8}~{}{\rm M}_{\odot})^{0.154}$
(Scott & Graham, 2013; Kormendy & Ho, 2013; Georgiev et al., 2016), and last
the inflow rate of an AGN disk in units of the Eddington rate varies within
${\dot{M}}_{\rm inflow}/{\dot{M}}_{\rm Edd}\sim 0.004-1.2$, where
${\dot{M}}_{\rm Edd}\equiv L_{\rm Edd}/c^{2}\eta_{\rm rad}$ with radiative
efficiency of $\eta_{\rm rad}=0.1$, and $L_{\rm Edd}$ is the Eddington
luminosity of the SMBH.

We find that the optical luminosity of the flare is lower than the kinetic
power of the jet ($f_{\rm op/kin}<1$) for all the pairs (Table 3). Considering
the bolometric correction of the breakout emission in the optical band ($\sim
10$, Paper I), our scenario demands $f_{\rm op/kin}\lesssim 0.1$. For the
pairs that do not satisfy this condition, we need to use $f_{\rm jet/BHL}>10$.
Note that a higher value ($f_{\rm jet/BHL}\lesssim 100$) is possible by
considering the enhancement of the accretion rate due to shocks caused by
recoil kicks at merger (Appendix. B of Paper I).

In the right panel of Fig. 2, the SEDs for thermal emission from the breakout
of the jet head for Pairs 1–12 are displayed. As they are bright at $\sim
10^{17}$–$10^{19}~{}{\rm Hz}$ and the duration is long ($\gtrsim 10^{6}~{}{\rm
s}$), they can be observed by X-ray telescopes, such as the Swift X-ray
telescope (XRT, Burrows et al. 2005), Chandra, XMM Newton (Jansen et al.,
2001), and the Nuclear Spectroscopic Telescope Array (NuSTAR, Harrison et al.
2013). Note that the shock cooling emission associated with the breakout
emission in all the pairs is so dim ($\ll 10^{40}~{}{\rm erg/s}$) that it is
buried within the AGN variability.

The derived values of $R_{\rm BH}/R_{\rm NSC}\sim 0.2-0.7\sim 1-8$ pc are much
larger than in the shock-cooling case. The different locations are driven by
the phases of the emission. In our model, the duration of flares corresponds
to the timescale at which the diffusion timescale is equal to the dynamical
timescale (e.g. Eqs. 10 and 2.3.1 for the shock cooling and breakout emission
cases, respectively). This timescale is inversely proportional to the density
of the ejecta or of the AGN disk, as well as to the square of the velocity of
the ejecta or of the shocks in the shock cooling and breakout emission
scenarios, respectively. As a result, since the velocities are typically
comparable (Tables 2 and 3), the density of the AGN disk in the breakout
emission scenario must be similar to the density of the ejecta in the shock
cooling emission scenario, in order to reproduce the duration of the flares.
Additionally, in the shock cooling emission scenario the density of the ejecta
is much lower than the local AGN density due to adiabatic expansion (Eq. 12).
Thus, the AGN density is typically much lower in the breakout emission
scenario, and therefore the distance from the SMBH needs to be larger. These
larger distance are consistent with the scenario that BHs in nuclear star
clusters are captured by AGN disks and merge with one other. Here, the
migration timescale of objects in AGN disks around massive SMBHs as inferred
for the flaring AGNs (Table 1) is so long that migration before BH mergers is
less expected.

To check whether the values derived for the aspect ratio ($h_{\rm AGN}$) are
plausible, we calculate the aspect ratio of the AGN disk ($h_{\rm AGN,TQM}$)
adopting the model in Thompson et al. (2005), using $R_{\rm BH}$, $M_{\rm
SMBH}$, ${\dot{M}}_{\rm inflow}$, and $m$. We find that the values of $h_{\rm
AGN,TQM}$ are roughly comparable to $h_{\rm AGN}$ within a factor of $\sim 2$
(9th column of Table 3). Such moderate differences could arise due to the
variation of $f_{\rm corr}$ reflecting the variation of the inclination of
jets. In Table 3, we list $f_{\rm corr}$ at which $h_{\rm AGN,TQM}=h_{\rm
AGN}$. This shows that $h_{\rm AGN,TQM}=h_{\rm AGN}$ is satisfied by
reasonable values for $f_{\rm corr}$ around $\sim 3$–$5$.

By comparing $L_{\rm AGN}/L_{\rm Edd}$ in Table 1 and ${\dot{M}}_{\rm
inflow}/{\dot{M}}_{\rm Edd}$ in Table 3, we can determine the value of $m$ at
which $L_{\rm AGN}/L_{\rm Edd}={\dot{M}}_{\rm inflow}/{\dot{M}}_{\rm Edd}$ as
${\dot{M}}_{\rm inflow}\propto m$. For the pairs 11 and 12, $m\geq 3$ is
required to satisfy $L_{\rm AGN}/L_{\rm Edd}\leq{\dot{M}}_{\rm
inflow}/{\dot{M}}_{\rm Edd}$. In the other pairs (1–10), $m$ can be as low as
$\geq 0.003$–$0.5$ to satisfy this condition. We presume that $m\lesssim 1$ is
acceptable (e.g. Hopkins & Quataert, 2010; Begelman & Armitage, 2023; Begelman
& Silk, 2023), although there are no studies constraining the possible ranges
of $m$ as far as we know. Note that the inflow rate at pc scales can be much
higher than the accretion rate onto the central SMBH (e.g. Izumi et al.,
2023), since a large fraction of gas is possibly consumed by star formation
and outflows, as predicted in Blandford & Begelman (1999), Thompson et al.
(2005), and Collin & Zahn (2008). Then ${\dot{M}}_{\rm inflow}/{\dot{M}}_{\rm
Edd}$ is allowed to be much larger than $L_{\rm AGN}/L_{\rm Edd}$, and the
required value of $m$ is enhanced. Hence, pairs 11 and 12 are not well
explained by breakout emission due to the required high values of $m$. On the
other hand, it is not clear whether AGN disks are usually in steady state,
which might complicate the comparison between the inflow rate and the AGN
luminosity.

##  4\. Tests of the Model

In the following we discuss how our model can be tested by examining the
distribution expected for the observable properties. From § 3.1 and 3.2 we see
that values for each model parameter are distributed in a narrow range (Tables
2 and 3). We discuss whether this is because flares originating from merging
BHs tend to have well-defined characteristic properties, or because of
observational selection effects in the way the search was conducted by Graham
et al. (2023). Flares were searched with specific ranges of parameters, and in
particular with the rise time within 5–100$~{}{\rm day}$, the decay time
within 10–200$~{}{\rm day}$, and the delay time from the GW detection in the
interval 0–200$~{}{\rm day}$.

In order for electromagnetic flares to be found in association with BH
mergers, BHs need to typically merge in bright AGNs. This is because most BH
mergers reported by LIGO/Virgo/KAGRA are found to merge at luminosity
distances of several Gpc (Abbott et al., 2020). At such large distances, AGNs
are easily missed unless they are bright. Indeed, the host SMBH masses for the
flares reported by Graham et al. (2023) are $\gtrsim 10^{8}~{}{\rm
M}_{\odot}$, hence so massive that AGNs are rarely missed in AGN searches at
the distances of the GW events. Also, assuming a luminosity distance of
$d_{\rm L}\sim 3~{}{\rm Gpc}$, SMBH mass of $10^{8}~{}{\rm M}_{\odot}$,
Eddington ratio of $\sim 1$, bolometric correction of 5, and the fraction of
the variable luminosity compared to the average luminosity $\sim 0.1$, the
flux of variable flares in the AGN is $\sim 2\times 10^{-13}~{}{\rm
erg/s/cm^{2}}$, which falls just above the sensitivity of ZTF $\sim
10^{-13}~{}{\rm erg/s/cm^{2}}$. Conversely, it is difficult to find flares
associated with BH mergers occurring in AGNs around less massive SMBHs through
surveys for AGN variability.

We next estimate the detection rate of electromagnetic flares associated with
BH mergers based on our scenario. If BH mergers actually produce
electromagnetic flares as found in Graham et al. (2023), the rate of such
flares is comparable to or less than the rate of BH mergers ($\sim
O(10)~{}{\rm Gpc^{-3}yr^{-1}}$). If flares within $\lesssim 3~{}{\rm Gpc}$ are
observable by the ZTF and all merging BHs produce electromagnetic flares, up
to $N_{\rm BBH,3Gpc}\sim 300$ flares associated with BH mergers can be
detected per year. On the other hand, we have estimated that a small fraction
of BH mergers ($f_{\rm EM/BBH}\lesssim 0.02$) could accompany detectable
electromagnetic flares in our model (see § 4.1 in Paper I for discussions).
Then, the seven flares discovered by (Graham et al., 2023) is roughly
consistent with the number of flares expected during LIGO/Virgo/KAGRA O3,
which is estimated as $N_{\rm EM/GW,O3}\sim N_{\rm BBH,3Gpc}f_{\rm
EM/BBH}t_{\rm O3}\lesssim 10$, where $t_{\rm O3}\sim 1.1$ is the duration of
the O3 operation in the unit of year. Constraints on the frequency of flares,
whose properties can be explained by emission from BHs, are hence useful to
test our model. Here, note that dimmer more frequent flares would be
contributed by solitary BHs (Tagawa et al., 2023c).

###  4.1. Shock cooling emission

From Table 2, we note that there are moderate variations in $R_{\rm BH}$ and
$T_{\rm BB}$ in the shock cooling emission scenario. To consider a possible
test of this model, we discuss its dependence on the observed properties, and
the range of these properties for which the shock cooling emission scenario is
inconsistent. If there are events with high or low $T_{\rm BB}$, extremely
large or small $R_{\rm BH}$, or $R_{\rm BO}/R_{\rm diff}\geq 1$, such flares
are inconsistent with the shock cooling emission scenario. However, the
dependence of the variables $T_{\rm BB}$, $R_{\rm BH}$, $R_{\rm BO}/R_{\rm
diff}$ on the ratio of the observed parameters ($L_{\rm AGN}$, $L_{\rm obs}$,
$m_{\rm BH}$, $M_{\rm SMBH}$, and $t_{\rm diff}$) over the adjustment
parameters ($f_{\rm corr}$ and $f_{\rm jet/BHL}$) is similar to the dependence
of $H_{\rm AGN}/R_{\rm BH}$ on the same ratio. For example, for $\beta_{\rm
c}/\theta_{0}<1$, $H_{\rm AGN}/R_{\rm BH}$ depends on the observed properties
as

$\displaystyle\frac{H_{\rm AGN}}{R_{\rm BH}}=\frac{9^{4/3}L_{\rm
AGN}^{2}L_{\rm obs}^{5}\kappa_{\rm ej}^{7}\theta_{0}^{2/3}f_{\rm
corr}^{9}}{2^{4/3}4^{4}\pi^{7}f_{\rm jet/BHL}^{4}f_{\rm
cons}^{2}c^{19}\alpha^{2}\eta_{\rm rad}^{2}A^{4/3}G^{3}}$
$\displaystyle\times\frac{1}{m_{\rm BH}^{8/3}M_{\rm SMBH}^{1/3}t_{\rm
diff}^{16/3}}$ $\displaystyle\propto L_{\rm obs}^{5}L_{\rm AGN}^{2}t_{\rm
delay}^{-16/3}m_{\rm BH}^{-8/3}f_{\rm corr}^{9}f_{\rm jet/BHL}^{-4},$ (39)

and $T_{\rm BB}$ depends as

$\displaystyle T_{\rm BB}=\frac{9^{3/8}L_{\rm AGN}^{3/8}L_{\rm
obs}^{11/8}\kappa_{\rm ej}^{3/2}f_{\rm
corr}^{9/4}}{(4\pi)^{1/3}(32\pi^{4})^{3/8}f_{\rm jet/BHL}^{9/8}f_{\rm
cons}^{3/8}c^{9/2}A^{3/8}G^{3/4}}$
$\displaystyle\times\frac{1}{\alpha^{3/8}\eta_{\rm rad}^{3/8}m_{\rm
BH}^{3/4}t_{\rm diff}^{5/4}a^{3/8}}$ $\displaystyle\propto L_{\rm
obs}^{11/8}L_{\rm AGN}^{3/8}t_{\rm delay}^{-5/4}m_{\rm BH}^{-3/4}f_{\rm
corr}^{9/4}f_{\rm jet/BHL}^{-9/8}.$ (40)

Then, if we adjust $f_{\rm corr}$ and $f_{\rm jet/BHL}$ so that $H_{\rm
AGN}/R_{BH}$ is consistent with a Shakura-Sunyaev disk, these parameters
($T_{\rm BB}$, $R_{\rm BH}$, and $R_{\rm BO}/R_{\rm diff}$) also fall in a
range of possible values (similar values to those derived in Table 2), even if
events with wide ranges of observed parameters ($L_{\rm AGN}$, $L_{\rm obs}$,
$m_{\rm BH}$, $M_{\rm SMBH}$, and $t_{\rm diff}$) are observed. Thus, in this
model, once $H_{\rm AGN}/R_{\rm BH}$ is adjusted to the value expected in the
Shakura-Sunyaev model, $T_{\rm BB}$, $R_{\rm BH}$, and $R_{\rm BO}/R_{\rm
diff}$ are then characterized by realistic values, and the distribution of the
observables becomes difficult to use as a further test of the model.

An interesting test of the model is the correlation between the delay time and
the duration of the flare. This is because the delay time is comparable to the
duration of a flare as long as $t_{\rm diff}\gg t_{\rm break}$, which is
satisfied in pairs 1–12 (Table 2). To constrain the correlation coefficient
e.g. with uncertainty of $\lesssim 0.3$ by 95 percentile, more than $\gtrsim
50$ events are needed to be observed. Hence, more events will be very helpful
as further diagnostics. Note that half of the pairs reported by Graham et al.
(2023) are unreal as a flare is associated with a few GW events, and hence, we
need to derive the correlation excluding the influence from false
associations. Also, we derived the delay time using digitizer. In addition,
the time dependence of the luminosity assumed by Graham et al. (2023) is
different from that expected in the shock cooling emission. Thus, both the
delay time and the duration suffer from significant uncertainties. If these
timescales can be well constrained, the model can be tested for each event by
comparing the delay time and the duration.

Another possible test is the detection of the cocoon breakout emission
preceding the shock cooling emission. In pairs 1–12, the duration and
temperature of the cocoon breakout emission range in the intervals $\sim
10$–$10^{5}~{}{\rm s}$ and $0.4$–$5~{}{\rm keV}$; thus wide-field X-ray
surveys, such as Einstein Probe (Yuan et al., 2015) and HiZ-GUNDAM (Yonetoku
et al., 2020) will be useful to detect the early breakout emission (left panel
of Fig. 2). Assuming that the luminosity of the breakout emission is similar
in magnitude to the jet power, the detectable distance is estimated to be

$\displaystyle d_{\rm det}=\left(\frac{L_{\rm jet}}{4\pi F_{\rm
sen}}\right)^{1/2}$ $\displaystyle\sim 1~{}{\rm Gpc}\left(\frac{L_{\rm
jet}}{10^{46}~{}{\rm erg/s}}\right)^{1/2}\left(\frac{F_{\rm
sen}}{10^{-11}~{}{\rm erg/s/cm^{2}}}\right)^{-1/2},$ (41)

where $F_{\rm sen}$ is the sensitivity of the facility. Here, the closest
distance among the flares discovered by Graham et al. (2023) is $\sim{\rm
Gpc}$. Since the duration of the cocoon breakout emission is $\sim
10$–$10^{5}~{}{\rm s}$ in pairs 1–12 and the sensitivity of Einstein Probe
(Yuan et al., 2015) and HiZ-GUNDAM is $F_{\rm sen}\sim 10^{-11}~{}{\rm
erg/s/cm^{2}}$ for $t_{\rm int}\sim 10^{4}~{}{\rm s}$, events at the
luminosity distance of $\sim{\rm Gpc}$ can be detected if $L_{\rm jet}\gtrsim
10^{46}~{}{\rm erg/s}$. If both the breakout emission and the shock cooling
emission are detected from the same AGN, it can be a robust test of this
model.

It is notable that in the shock cooling model, the color keeps evolving to
redder in all pairs (Fig. 3 d). The $m_{g}-m_{r}$ evolves by $\sim 0.4$ in
$50~{}{\rm days}$. Such mild evolution of the color can be a strong test of
this model. On the other hand, for non-thermal emission from the breakout of
the jet head, the color is expected to be unchanged. Note that the colors are
presented for pair 3 in Fig. 2 of Graham et al. (2020), revealing almost no
evolution. However, due to the contamination from the host AGN emission, any
color evolution is likely difficult to constrain. To derive the color
evolution more precisely, observations at additional frequencies, e.g. by
ULTRASAT (black lines of Fig. 3 c), would also be useful. Hence, when
available, the evolution of the color will be an additional diagnostic in
future observations.

###  4.2. Breakout emission

Next we discuss ways to test the breakout emission scenario. To do this, we
first present the dependence of physical quantities on observables as

$\displaystyle\beta_{\rm h}\propto\theta_{0}^{-12/11}f_{\rm
jet/BHL}^{3/11}f_{\rm corr}^{-3/11}t_{\rm delay}^{-3/11}t_{\rm
duration}^{1/22},$ (42) $\displaystyle H_{\rm
AGN}\propto\theta_{0}^{-12/11}f_{\rm jet/BHL}^{3/11}f_{\rm
corr}^{-14/11}t_{\rm delay}^{8/11}t_{\rm duration}^{1/22},$ (43)
$\displaystyle\rho_{\rm AGN}\propto\theta_{0}^{24/11}f_{\rm
jet/BHL}^{-6/11}f_{\rm corr}^{6/11}t_{\rm delay}^{6/11}t_{\rm
duration}^{-12/22},$ (44) $\displaystyle R_{\rm
BH}\propto\theta_{0}^{-8/11}f_{\rm jet/BHL}^{2/11}M_{\rm SMBH}^{1/3}f_{\rm
corr}^{-2/11}t_{\rm delay}^{-2/11}t_{\rm duration}^{4/22},$ (45)
$\displaystyle L_{\rm j}\propto\theta_{0}^{8/11}f_{\rm jet/BHL}^{9/11}f_{\rm
corr}^{-9/11}t_{\rm delay}^{13/11}t_{\rm duration}^{-19/22},$ (46)
$\displaystyle{\dot{M}}_{\rm inflow}\propto\theta_{0}^{4/11}f_{\rm
jet/BHL}^{-1/11}mM_{\rm SMBH}^{1/2}f_{\rm corr}^{-17/11}t_{\rm
delay}^{23/11}t_{\rm duration}^{-13/22}.$ (47)

For the breakout emission from the jet head, if the shocked gas becomes
relativistic ($\beta_{\rm h}\gtrsim 0.8$), the probability of observing non-
thermal emission is significantly reduced due to the relativistic beaming
effects, the shift of the minimum energy, and the time dilution (as shown in
Fig. 3 b of Paper I). Thus, if the shock is relativistic, the breakout
emission is presumably not observed by current facilities, such as ZTF. Since
$\beta_{\rm h}$ is estimated to be $\sim 0.2$–$0.4$ for the observed flares,
if $\beta_{\rm h}$ is higher by a factor of $\sim 3$ compared to our
estimates, flares cannot be explained by breakout emission. $\beta_{\rm
h}\gtrsim 0.8$ is satisfied in all the pairs if events with delay time of
$t_{\rm delay}\lesssim 2~{}{\rm day}$ are found. Note that this is not
compensated by $f_{\rm corr}$ and $f_{\rm jet/BHL}$ (as in the shock cooling
emission scenario). This is because $\beta_{\rm h}$ is reduced only by a
factor of $\sim 1.1$ by enhancing $f_{\rm corr}$ to the maximum value ($\sim
1/\theta_{0}$). Also, if $f_{\rm jet/BHL}$ is reduced to lower $\beta_{\rm
h}$, $f_{\rm op/kin}$ becomes larger than 1, which violates another
requirement for the model. If associations between the flares and GWs are due
to random coincidence, the delay time is expected to be distributed uniformly
in the range of $0$–$200~{}{\rm day}$. Assuming a uniform distribution, we can
test the breakout emission scenario at the $\sim 1~{}\sigma$ and $\sim
2~{}\sigma$ levels after discovering $\gtrsim 100$ and $\gtrsim 300$ events,
respectively, by checking if there are events with $t_{\rm delay}\lesssim
2~{}{\rm day}$. If we find optical flares with $t_{\rm delay}\lesssim 2~{}{\rm
day}$, the breakout emission scenario is disfavored.

If $t_{\rm duration}$ is enhanced by one order of magnitude, the jet power is
reduced by a similar factor, and the fraction of the optical luminosity to the
jet power exceeds one for Pairs 1–3, 6, 11, and 12, and then, the breakout
emission scenario becomes inconsistent for these pairs. However, the
enhancement of $L_{j}$ due to long $t_{\rm duration}$ can be compensated by
$f_{\rm jet}$ and $f_{\rm corr}$ up to about a factor of $\sim 10$, which is
limited by the requirement for $\beta_{\rm h}$ $(\lesssim 0.8)$ as discussed
above. Hence, the model is not well tested by $t_{\rm duration}$.

Currently, the properties of the observed flares satisfy the conditions
required for the breakout emission scenario. To be consistent with this
scenario, the delay time should not be shorter than $\sim 2~{}{\rm day}$, and
the color of flares should be distributed in a narrow range as discussed in §
3.1. Thus, to test whether there are flares with properties being inconsistent
with the breakout emission scenario, more events will need to be observed.

##  5\. Conclusions

In this paper we have presented the properties of emission from shocks
emerging from collisions between AGN gas and a jet launched from a merger
remnant BH in an AGN disk. Our model includes the evolution throughout the
shock breakout and subsequent cooling emission phases. We then applied this
model to the candidate flares reported in Graham et al. (2023). Our results
are summarized as follows.

  1. 1.

We fit the characteristic features of all of the events with both emission
processes. Both processes could fit each observation, suggesting that such
fits may themselves be insufficient to rule out the AGN origin for the flares.
The reconstructed parameters might then be indicative of the selection
algorithm determining the false alarm rate of associations.

  2. 2.

While both processes could be made consistent with the observed events with
appropriate parameter selection, we found that the implied merger distance
from the central SMBH is markedly different for the two processes.
Specifically, shock cooling emission can explain the observed properties if
the mergers happen $R_{\rm BH}\sim 0.06$–$0.3~{}{\rm pc}$ from the SMBH, which
is consistent with the locations of AGN-assisted merger models (e.g. Tagawa et
al. 2020b). On the other hand, breakout emission would require a much larger
distance of $R_{\rm BH}\sim 1$–$8~{}{\rm pc}$ to explain the observed flare
duration. This may be possible for mergers in AGNs with high SMBH masses, in
which migration is inefficient (e.g. Perna et al. 2021b).

  3. 3.

Follow-up observations could help further constrain the reconstructed
parameters of the events. X-ray observations would have to be made prior to
the optical detection, which likely requires future wide-field surveys. In
addition, follow-up observations determining the spectral evolution of the
electromagnetic flares would be important.

H.T. was supported by the National Key R&D Program of China (Grant No.
2021YFC2203002) and the National Natural Science Foundation of China (Grant
No. 12173071). S.S.K. was supported by Japan Society for the Promotion of
Science (JSPS) KAKENHI grant Number 22K14028, 21H04487, and 23H04899, and the
Tohoku Initiative for Fostering Global Researchers for Interdisciplinary
Sciences (TI-FRIS) of MEXT’s Strategic Professional Development Program for
Young Researchers. Z.H. was supported by NASA grant 80NSSC22K0822 and NSF
grant AST-2006176. R.P. acknowledges support by NSF award AST-2006839. I.B.
acknowledges the support of the Alfred P. Sloan Foundation and NSF grants
PHY-1911796 and PHY-2110060.

## References

  * Aasi et al. (2015) Aasi, J., et al. 2015, CQGra, 32, 074001, doi: 10.1088/0264-9381/32/7/074001
  * Abbott et al. (2020) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, arXiv e-prints, arXiv:2010.14527.  https://arxiv.org/abs/2010.14527
  * Acernese et al. (2015) Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, Classical and Quantum Gravity, 32, 024001, doi: 10.1088/0264-9381/32/2/024001
  * Akutsu et al. (2021) Akutsu, T., Ando, M., Arai, K., et al. 2021, Progress of Theoretical and Experimental Physics, 2021, 05A101, doi: 10.1093/ptep/ptaa125
  * Antonini et al. (2022) Antonini, F., Gieles, M., Dosopoulou, F., & Chattopadhyay, D. 2022, arXiv e-prints, arXiv:2208.01081, doi: 10.48550/arXiv.2208.01081
  * Antonini et al. (2017) Antonini, F., Toonen, S., & Hamers, A. S. 2017, ApJ, 841, 77, doi: 10.3847/1538-4357/aa6f5e
  * Arca Sedda et al. (2023) Arca Sedda, M., Naoz, S., & Kocsis, B. 2023, arXiv e-prints, arXiv:2302.14071, doi: 10.48550/arXiv.2302.14071
  * Arnett (1980) Arnett, W. D. 1980, ApJ, 237, 541, doi: 10.1086/157898
  * Ashton et al. (2020) Ashton, G., Ackley, K., Magaña Hernandez, I., & Piotrzkowski, B. 2020, arXiv e-prints, arXiv:2009.12346.  https://arxiv.org/abs/2009.12346
  * Bagoly et al. (2016) Bagoly, Z., Szécsi, D., Balázs, L. G., et al. 2016, A&A, 593, L10, doi: 10.1051/0004-6361/201628569
  * Banerjee (2017) Banerjee, S. 2017, MNRAS, 467, 524.  http://adsabs.harvard.edu/abs/2017MNRAS.467..524B
  * Barack et al. (2019) Barack, L., Cardoso, V., Nissanke, S., et al. 2019, Classical and Quantum Gravity, 36, 143001, doi: 10.1088/1361-6382/ab0587
  * Bartos et al. (2017) Bartos, I., Kocsis, B., Haiman, Z., & Márka, S. 2017, ApJ, 835, 165\.  http://adsabs.harvard.edu/abs/2017ApJ...835..165B
  * Bartos et al. (2023) Bartos, I., Rosswog, S., Gayathri, V., et al. 2023, arXiv e-prints, arXiv:2302.10350, doi: 10.48550/arXiv.2302.10350
  * Begelman & Armitage (2023) Begelman, M. C., & Armitage, P. J. 2023, MNRAS, 521, 5952, doi: 10.1093/mnras/stad914
  * Begelman & Silk (2023) Begelman, M. C., & Silk, J. 2023, arXiv e-prints, arXiv:2305.19081, doi: 10.48550/arXiv.2305.19081
  * Belczynski et al. (2016) Belczynski, K., Daniel, E. H., Bulik, T., & O’Shaughnessy, R. 2016, Nature, 534, 512.  http://adsabs.harvard.edu/abs/2016Natur.534..512B
  * Bellovary et al. (2016) Bellovary, J. M., Mac Low, M. M., McKernan, B., & Ford, K. E. S. 2016, ApJ, 819, L17.  http://adsabs.harvard.edu/abs/2016ApJ...819L..17B
  * Berger (2014) Berger, E. 2014, ARA&A, 52, 43, doi: 10.1146/annurev-astro-081913-035926
  * Blandford & Begelman (1999) Blandford, R. D., & Begelman, M. C. 1999, MNRAS, 303, L1, doi: 10.1046/j.1365-8711.1999.02358.x
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433, doi: 10.1093/mnras/179.3.433
  * Boekholt et al. (2023) Boekholt, T. C. N., Rowan, C., & Kocsis, B. 2023, MNRAS, 518, 5653, doi: 10.1093/mnras/stac3495
  * Bromberg et al. (2011) Bromberg, O., Nakar, E., Piran, T., & Sari, R. 2011, ApJ, 740, 100, doi: 10.1088/0004-637X/740/2/100
  * Burrows et al. (2005) Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, Space Sci. Rev., 120, 165, doi: 10.1007/s11214-005-5097-2
  * Calderón Bustillo et al. (2021) Calderón Bustillo, J., Leong, S. H. W., Chandra, K., McKernan, B., & Ford, K. E. S. 2021, arXiv e-prints, arXiv:2112.12481.  https://arxiv.org/abs/2112.12481
  * Chen et al. (2023a) Chen, K., Ren, J., & Dai, Z.-G. 2023a, ApJ, 948, 136, doi: 10.3847/1538-4357/acc45f
  * Chen et al. (2023b) Chen, Y.-X., Jiang, Y.-F., Goodman, J., & Ostriker, E. C. 2023b, ApJ, 948, 120, doi: 10.3847/1538-4357/acc023
  * Choksi et al. (2023) Choksi, N., Chiang, E., Fung, J., & Zhu, Z. 2023, MNRAS, doi: 10.1093/mnras/stad2269
  * Collin & Zahn (2008) Collin, S., & Zahn, J. P. 2008, A&A, 477, 419, doi: 10.1051/0004-6361:20078191
  * Connaughton et al. (2016) Connaughton, V., Burns, E., Goldstein, A., et al. 2016, ApJ, 826, L6, doi: 10.3847/2041-8205/826/1/L6
  * de Mink & King (2017) de Mink, S. E., & King, A. 2017, ApJ, 839, L7, doi: 10.3847/2041-8213/aa67f3
  * DeLaurentiis et al. (2023) DeLaurentiis, S., Haiman, Z., & Epstein-Martin, M. 2023, in American Astronomical Society Meeting Abstracts, Vol. 55, American Astronomical Society Meeting Abstracts, 351.02
  * Derdzinski et al. (2021) Derdzinski, A., D’Orazio, D., Duffell, P., Haiman, Z., & MacFadyen, A. 2021, MNRAS, 501, 3540, doi: 10.1093/mnras/staa3976
  * Derdzinski & Mayer (2023) Derdzinski, A., & Mayer, L. 2023, MNRAS, 521, 4522, doi: 10.1093/mnras/stad749
  * Derdzinski et al. (2019) Derdzinski, A. M., D’Orazio, D., Duffell, P., Haiman, Z., & MacFadyen, A. 2019, MNRAS, 486, 2754, doi: 10.1093/mnras/stz1026
  * Dominik et al. (2012) Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52.  http://adsabs.harvard.edu/abs/2012ApJ...759...52D
  * Duras et al. (2020) Duras, F., Bongiorno, A., Ricci, F., et al. 2020, A&A, 636, A73, doi: 10.1051/0004-6361/201936817
  * Ford & McKernan (2022) Ford, K. E. S., & McKernan, B. 2022, MNRAS, 517, 5827, doi: 10.1093/mnras/stac2861
  * Fragione & Kocsis (2018) Fragione, G., & Kocsis, B. 2018, Phys. Rev. Lett., 121, 161103, doi: 10.1103/PhysRevLett.121.161103
  * Fragione & Kocsis (2019) —. 2019, MNRAS, 486, 4781, doi: 10.1093/mnras/stz1175
  * Gayathri et al. (2023) Gayathri, V., Wysocki, D., Yang, Y., et al. 2023, arXiv e-prints, arXiv:2301.04187, doi: 10.48550/arXiv.2301.04187
  * Gayathri et al. (2021) Gayathri, V., Yang, Y., Tagawa, H., Haiman, Z., & Bartos, I. 2021, ApJ, 920, L42, doi: 10.3847/2041-8213/ac2cc1
  * Gayathri et al. (2022) Gayathri, V., Healy, J., Lange, J., et al. 2022, Nature Astronomy, 6, 344, doi: 10.1038/s41550-021-01568-w
  * Generozov & Perets (2022) Generozov, A., & Perets, H. B. 2022, arXiv e-prints, arXiv:2212.11301, doi: 10.48550/arXiv.2212.11301
  * Georgiev et al. (2016) Georgiev, I. Y., Boker, T., Leigh, N., Lutzgendorf, N., & Neumayer, N. 2016, MNRAS, 457, 2122.  http://adsabs.harvard.edu/abs/2016MNRAS.457.2122G
  * Goodman & Tan (2004) Goodman, J., & Tan, J. C. 2004, ApJ, 608, 108.  http://adsabs.harvard.edu/abs/2004ApJ...608..108G
  * Gottlieb et al. (2018) Gottlieb, O., Nakar, E., & Piran, T. 2018, MNRAS, 473, 576, doi: 10.1093/mnras/stx2357
  * Graham et al. (2020) Graham, M. J., Ford, K. E. S., McKernan, B., et al. 2020, Phys. Rev. Lett., 124, 251102, doi: 10.1103/PhysRevLett.124.251102
  * Graham et al. (2023) Graham, M. J., McKernan, B., Ford, K. E. S., et al. 2023, ApJ, 942, 99, doi: 10.3847/1538-4357/aca480
  * Grishin et al. (2021) Grishin, E., Bobrick, A., Hirai, R., Mandel, I., & Perets, H. B. 2021, MNRAS, 507, 156, doi: 10.1093/mnras/stab1957
  * Grishin et al. (2023) Grishin, E., Gilbaum, S., & Stone, N. C. 2023, arXiv e-prints, arXiv:2307.07546, doi: 10.48550/arXiv.2307.07546
  * Hada et al. (2013) Hada, K., Kino, M., Doi, A., et al. 2013, ApJ, 775, 70, doi: 10.1088/0004-637X/775/1/70
  * Hada et al. (2018) Hada, K., Doi, A., Wajima, K., et al. 2018, ApJ, 860, 141, doi: 10.3847/1538-4357/aac49f
  * Haiman et al. (2009) Haiman, Z., Kocsis, B., & Menou, K. 2009, ApJ, 700, 1952, doi: 10.1088/0004-637X/700/2/1952
  * Hamidani & Ioka (2023) Hamidani, H., & Ioka, K. 2023, MNRAS, 524, 4841, doi: 10.1093/mnras/stad1933
  * Harrison et al. (2013) Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103, doi: 10.1088/0004-637X/770/2/103
  * Hopkins & Quataert (2010) Hopkins, P. F., & Quataert, E. 2010, MNRAS, 407, 1529, doi: 10.1111/j.1365-2966.2010.17064.x
  * Izumi et al. (2023) Izumi, T., Wada, K., Imanishi, M., et al. 2023, arXiv e-prints, arXiv:2305.03993, doi: 10.48550/arXiv.2305.03993
  * Jansen et al. (2001) Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1, doi: 10.1051/0004-6361:20000036
  * Kimura et al. (2021) Kimura, S. S., Murase, K., & Bartos, I. 2021, ApJ, 916, 111, doi: 10.3847/1538-4357/ac0535
  * Kinugawa et al. (2014) Kinugawa, T., Inayoshi, K., Hotokezaka, K., Nakauchi, D., & T., N. 2014, MNRAS, 442, 2963.  http://adsabs.harvard.edu/abs/2014MNRAS.442.2963K
  * Kitaki et al. (2021) Kitaki, T., Mineshige, S., Ohsuga, K., & Kawashima, T. 2021, PASJ, 73, 450, doi: 10.1093/pasj/psab011
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Kumamoto et al. (2018) Kumamoto, J., Fujii, M. S., & Tanikawa, A. 2018, arXiv e-prints.  https://arxiv.org/abs/1811.06726
  * Lazzati et al. (2022) Lazzati, D., Soares, G., & Perna, R. 2022, ApJ, 938, L18, doi: 10.3847/2041-8213/ac98ad
  * Levin & Beloborodov (2003) Levin, Y., & Beloborodov, A. M. 2003, ApJ, 590, L33, doi: 10.1086/376675
  * Li et al. (2023) Li, J., Dempsey, A. M., Li, H., Lai, D., & Li, S. 2023, ApJ, 944, L42, doi: 10.3847/2041-8213/acb934
  * Li et al. (2022) Li, J., Lai, D., & Rodet, L. 2022, ApJ, 934, 154, doi: 10.3847/1538-4357/ac7c0d
  * Li et al. (2021) Li, Y.-P., Dempsey, A. M., Li, S., Li, H., & Li, J. 2021, ApJ, 911, 124, doi: 10.3847/1538-4357/abed48
  * Liska et al. (2018) Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2018, MNRAS, 474, L81, doi: 10.1093/mnrasl/slx174
  * Liska et al. (2021) —. 2021, MNRAS, 507, 983, doi: 10.1093/mnras/staa099
  * Martinez et al. (2022) Martinez, M. A. S., Rodriguez, C. L., & Fragione, G. 2022, ApJ, 937, 78, doi: 10.3847/1538-4357/ac8d55
  * McKernan et al. (2018) McKernan, B., Ford, K. E. S., Bellovary, J., et al. 2018, ApJ, 866, 66, doi: 10.3847/1538-4357/aadae5
  * McKernan et al. (2019) McKernan, B., Ford, K. E. S., Bartos, I., et al. 2019, ApJ, 884, L50, doi: 10.3847/2041-8213/ab4886
  * Michaely & Perets (2019) Michaely, E., & Perets, H. B. 2019, ApJ, 887, L36, doi: 10.3847/2041-8213/ab5b9b
  * Milosavljević & Loeb (2004) Milosavljević, M., & Loeb, A. 2004, ApJ, 604, L45, doi: 10.1086/383467
  * Miralda-Escudé & Kollmeier (2005) Miralda-Escudé, J., & Kollmeier, J. A. 2005, ApJ, 619, 30, doi: 10.1086/426467
  * Morag et al. (2023) Morag, J., Sapir, N., & Waxman, E. 2023, MNRAS, 522, 2764, doi: 10.1093/mnras/stad899
  * Morton et al. (2023) Morton, S., Rinaldi, S., Torres-Orjuela, A., et al. 2023, arXiv e-prints, arXiv:2310.16025.  https://arxiv.org/abs/2310.16025
  * Nakar & Piran (2017) Nakar, E., & Piran, T. 2017, ApJ, 834, 28, doi: 10.3847/1538-4357/834/1/28
  * Nakar & Sari (2010) Nakar, E., & Sari, R. 2010, ApJ, 725, 904, doi: 10.1088/0004-637X/725/1/904
  * Narayan et al. (2021) Narayan, R., Chael, A., Chatterjee, K., Ricarte, A., & Curd, B. 2021, arXiv e-prints, arXiv:2108.12380.  https://arxiv.org/abs/2108.12380
  * O’Leary et al. (2016) O’Leary, R. M., Meiron, Y., & Kocsis, B. 2016, ApJ, 824, L12, doi: 10.3847/2041-8205/824/1/L12
  * Ostriker (1983) Ostriker, J. P. 1983, ApJ, 273, 99, doi: 10.1086/161351
  * Palmese et al. (2021) Palmese, A., Fishbach, M., Burke, C. J., Annis, J., & Liu, X. 2021, ApJ, 914, L34, doi: 10.3847/2041-8213/ac0883
  * Panaitescu & Kumar (2001) Panaitescu, A., & Kumar, P. 2001, ApJ, 560, L49, doi: 10.1086/324061
  * Perna et al. (2021a) Perna, R., Lazzati, D., & Cantiello, M. 2021a, ApJ, 906, L7, doi: 10.3847/2041-8213/abd319
  * Perna et al. (2021b) Perna, R., Tagawa, H., Haiman, Z., & Bartos, I. 2021b, ApJ, 915, 10, doi: 10.3847/1538-4357/abfdb4
  * Perna et al. (2019) Perna, R., Wang, Y.-H., Farr, W. M., Leigh, N., & Cantiello, M. 2019, ApJ, 878, L1, doi: 10.3847/2041-8213/ab2336
  * Piro et al. (2021) Piro, A. L., Haynie, A., & Yao, Y. 2021, ApJ, 909, 209, doi: 10.3847/1538-4357/abe2b1
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Polko & McKinney (2017) Polko, P., & McKinney, J. C. 2017, MNRAS, 464, 2660, doi: 10.1093/mnras/stw1875
  * Portegies Zwart & McMillan (2000) Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJ, 528, L17.  http://adsabs.harvard.edu/abs/2000ApJ...528L..17P
  * Prasad et al. (2023) Prasad, C., Wang, Y., Perna, R., Ford, K. E. S., & McKernan, B. 2023, arXiv e-prints, arXiv:2310.00020, doi: 10.48550/arXiv.2310.00020
  * Qian et al. (2023) Qian, K., Li, J., & Lai, D. 2023, arXiv e-prints, arXiv:2310.12208, doi: 10.48550/arXiv.2310.12208
  * Rasskazov & Kocsis (2019) Rasskazov, A., & Kocsis, B. 2019, ApJ, 881, 20, doi: 10.3847/1538-4357/ab2c74
  * Ray et al. (2023) Ray, M., Lazzati, D., & Perna, R. 2023, MNRAS, 521, 4233, doi: 10.1093/mnras/stad816
  * Rodriguez et al. (2016) Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016, Phys. Rev. D., 93, 084029.  http://adsabs.harvard.edu/abs/2016PhRvD..93h4029R
  * Rodríguez-Ramírez et al. (2023) Rodríguez-Ramírez, J. C., Bom, C. R., Fraga, B., & Nemmen, R. 2023, arXiv e-prints, arXiv:2304.10567, doi: 10.48550/arXiv.2304.10567
  * Romero-Shaw et al. (2022) Romero-Shaw, I., Lasky, P. D., & Thrane, E. 2022, ApJ, 940, 171, doi: 10.3847/1538-4357/ac9798
  * Romero-Shaw et al. (2023) Romero-Shaw, I. M., Gerosa, D., & Loutrel, N. 2023, MNRAS, 519, 5352, doi: 10.1093/mnras/stad031
  * Romero-Shaw et al. (2020) Romero-Shaw, I. M., Lasky, P. D., Thrane, E., & Calderon Bustillo, J. 2020, arXiv e-prints, arXiv:2009.04771.  https://arxiv.org/abs/2009.04771
  * Rossi et al. (2010) Rossi, E. M., Lodato, G., Armitage, P. J., Pringle, J. E., & King, A. R. 2010, MNRAS, 401, 2021, doi: 10.1111/j.1365-2966.2009.15802.x
  * Rowan et al. (2022) Rowan, C., Boekholt, T., Kocsis, B., & Haiman, Z. 2022, arXiv e-prints, arXiv:2212.06133, doi: 10.48550/arXiv.2212.06133
  * Rowan et al. (2023) Rowan, C., Whitehead, H., Boekholt, T., Kocsis, B., & Haiman, Z. 2023, arXiv e-prints, arXiv:2309.14433, doi: 10.48550/arXiv.2309.14433
  * Rozner et al. (2023) Rozner, M., Generozov, A., & Perets, H. B. 2023, MNRAS, doi: 10.1093/mnras/stad603
  * Samsing et al. (2014) Samsing, J., MacLeod, M., & Ramirez-Ruiz, E. 2014, ApJ, 784, 71, doi: 10.1088/0004-637X/784/1/71
  * Samsing et al. (2020) Samsing, J., Bartos, I., D’Orazio, D. J., et al. 2020, arXiv e-prints, arXiv:2010.09765.  https://arxiv.org/abs/2010.09765
  * Sapir & Waxman (2017) Sapir, N., & Waxman, E. 2017, ApJ, 838, 130, doi: 10.3847/1538-4357/aa64df
  * Scott & Graham (2013) Scott, N., & Graham, A. W. 2013, ApJ, 763, 76.  http://adsabs.harvard.edu/abs/2013ApJ...763...76S
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A& A, 24, 337.  http://adsabs.harvard.edu/abs/1973A%26A....24..337S
  * Shvartzvald et al. (2023) Shvartzvald, Y., Waxman, E., Gal-Yam, A., et al. 2023, arXiv e-prints, arXiv:2304.14482, doi: 10.48550/arXiv.2304.14482
  * Silsbee & Tremaine (2017) Silsbee, K., & Tremaine, S. 2017, ApJ, 836, 39.  http://adsabs.harvard.edu/abs/2017ApJ...836...39S
  * Spera et al. (2019) Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889, doi: 10.1093/mnras/stz359
  * Stone et al. (2017) Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, MNRAS, 464, 946.  http://adsabs.harvard.edu/abs/2017MNRAS.464..946S
  * Tagawa et al. (2020a) Tagawa, H., Haiman, Z., Bartos, I., & Kocsis, B. 2020a, ApJ, 899, 26, doi: 10.3847/1538-4357/aba2cc
  * Tagawa et al. (2020b) Tagawa, H., Haiman, Z., & Kocsis, B. 2020b, ApJ, 898, 25, doi: 10.3847/1538-4357/ab9b8c
  * Tagawa et al. (2023a) Tagawa, H., Kimura, S. S., & Haiman, Z. 2023a, ApJ, 955, 23, doi: 10.3847/1538-4357/ace71d
  * Tagawa et al. (2023b) Tagawa, H., Kimura, S. S., Haiman, Z., Perna, R., & Bartos, I. 2023b, ApJ, 950, 13, doi: 10.3847/1538-4357/acc4bb
  * Tagawa et al. (2023c) —. 2023c, ApJ, 946, L3, doi: 10.3847/2041-8213/acc103
  * Tagawa et al. (2022) Tagawa, H., Kimura, S. S., Haiman, Z., et al. 2022, ApJ, 927, 41, doi: 10.3847/1538-4357/ac45f8
  * Tagawa et al. (2021a) Tagawa, H., Kocsis, B., Haiman, Z., et al. 2021a, ApJ, 908, 194, doi: 10.3847/1538-4357/abd555
  * Tagawa et al. (2021b) —. 2021b, ApJ, 907, L20, doi: 10.3847/2041-8213/abd4d3
  * Tagawa et al. (2018) Tagawa, H., Kocsis, B., & Saitoh, R. T. 2018, Phys. Rev. Lett., 120, 261101\.  http://adsabs.harvard.edu/abs/2018PhRvL.120z1101T
  * Tanigawa & Tanaka (2016) Tanigawa, T., & Tanaka, H. 2016, ApJ, 823, 48, doi: 10.3847/0004-637X/823/1/48
  * Tanigawa & Watanabe (2002) Tanigawa, T., & Watanabe, S.-i. 2002, ApJ, 580, 506, doi: 10.1086/343069
  * Tanikawa et al. (2022) Tanikawa, A., Yoshida, T., Kinugawa, T., et al. 2022, ApJ, 926, 83, doi: 10.3847/1538-4357/ac4247
  * Tchekhovskoy et al. (2010) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2010, ApJ, 711, 50, doi: 10.1088/0004-637X/711/1/50
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. 2021, arXiv e-prints, arXiv:2111.03606.  https://arxiv.org/abs/2111.03606
  * Thompson et al. (2005) Thompson, T. A., Quataert, E., & Murray, N. 2005, ApJ, 630, 167.  http://adsabs.harvard.edu/abs/2005ApJ...630..167T
  * Wang et al. (2021a) Wang, J.-M., Liu, J.-R., Ho, L. C., & Du, P. 2021a, ApJ, 911, L14, doi: 10.3847/2041-8213/abee81
  * Wang et al. (2021b) Wang, J.-M., Liu, J.-R., Ho, L. C., Li, Y.-R., & Du, P. 2021b, ApJ, 916, L17, doi: 10.3847/2041-8213/ac0b46
  * Wang et al. (2023) Wang, Y., Zhu, Z., & Lin, D. N. C. 2023, arXiv e-prints, arXiv:2308.09129, doi: 10.48550/arXiv.2308.09129
  * Wang et al. (2022) Wang, Y.-H., Lazzati, D., & Perna, R. 2022, MNRAS, 516, 5935, doi: 10.1093/mnras/stac1968
  * Whitehead et al. (2023) Whitehead, H., Rowan, C., Boekholt, T., & Kocsis, B. 2023, arXiv e-prints, arXiv:2309.11561, doi: 10.48550/arXiv.2309.11561
  * Xin et al. (2023) Xin, C., Haiman, Z., Perna, R., Wang, Y., & Ryu, T. 2023, arXiv e-prints, arXiv:2303.12846, doi: 10.48550/arXiv.2303.12846
  * Yang et al. (2021) Yang, Y., Bartos, I., Fragione, G., et al. 2021, arXiv e-prints, arXiv:2105.02342.  https://arxiv.org/abs/2105.02342
  * Yonetoku et al. (2020) Yonetoku, D., Mihara, T., Doi, A., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11444, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114442Z, doi: 10.1117/12.2560603
  * Yuan et al. (2021) Yuan, C., Murase, K., Guetta, D., et al. 2021, arXiv e-prints, arXiv:2112.07653.  https://arxiv.org/abs/2112.07653
  * Yuan & Narayan (2014) Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529, doi: 10.1146/annurev-astro-082812-141003
  * Yuan et al. (2015) Yuan, W., Zhang, C., Feng, H., et al. 2015, arXiv e-prints, arXiv:1506.07735.  https://arxiv.org/abs/1506.07735
  * Zhou & Wang (2023) Zhou, Z.-H., & Wang, K. 2023, arXiv e-prints, arXiv:2310.15832.  https://arxiv.org/abs/2310.15832
  * Zhu (2023) Zhu, J.-P. 2023, arXiv e-prints, arXiv:2310.14255.  https://arxiv.org/abs/2310.14255
  * Zhu et al. (2021a) Zhu, J.-P., Wang, K., & Zhang, B. 2021a, ApJ, 917, L28, doi: 10.3847/2041-8213/ac1a17
  * Zhu et al. (2021b) Zhu, J.-P., Yang, Y.-P., Zhang, B., et al. 2021b, ApJ, 914, L19, doi: 10.3847/2041-8213/abff5a
  * Zhu et al. (2021c) Zhu, J.-P., Zhang, B., Yu, Y.-W., & Gao, H. 2021c, ApJ, 906, L11, doi: 10.3847/2041-8213/abd412