# Observable signature of merging stellar-mass black holes in active galactic
nuclei

Hiromichi Tagawa11affiliationmark:  22affiliationmark:  , Shigeo S
Kimura22affiliationmark:  33affiliationmark:  , Zoltán
Haiman11affiliationmark:  , Rosalba Perna44affiliationmark:
55affiliationmark:  , Imre Bartos66affiliationmark:  11affiliationmark:
Department of Astronomy, Columbia University, 550 W. 120th St., New York, NY,
10027, USA.
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 Physics and Astronomy, Stony Brook
University, Stony Brook, NY 11794-3800, USA.
55affiliationmark:  Center for Computational Astrophysics, Flatiron Institute,
New York, NY 10010, USA.
66affiliationmark:  Department of Physics, University of Florida, PO Box
118440, Gainesville, FL 32611, USA.  E-mail<EMAIL_ADDRESS>
###### Abstract

The origin of stellar-mass black hole mergers discovered through gravitational
waves is being widely debated. Mergers in the disks of active galactic nuclei
(AGN) represent a promising source of origin, with possible observational
clues in the gravitational wave data. Beyond gravitational waves, a unique
signature of AGN-assisted mergers is electromagnetic emission from the
accreting black holes. Here we show that jets launched by accreting black
holes merging in an AGN disk can be detected as peculiar transients by
infrared, optical, and X-ray observatories. We further show that this emission
mechanism can explain the possible associations between gravitational wave
events and the optical transient ZTF19abanrhr and the proposed gamma-ray
counterparts GW150914-GBM and LVT151012-GBM. We demonstrate how these
associations, if genuine, can be used to reconstruct the properties of these
events’ environments. Searching for infrared and X-ray counterparts to similar
electromagnetic transients in the future, once host galaxies are localized by
optical observations, could provide a smoking gun signature of the mergers’
AGN origin.

##  1\. Introduction

Despite the large number of black hole mergers discovered by the LIGO (LIGO
Scientific Collaboration et al., 2015), Virgo (Acernese et al., 2015) and
KAGRA (Akutsu et al., 2021) gravitational wave observatories (The LIGO
Scientific Collaboration et al., 2021), the astrophysical pathways to these
mergers remain debated. A promising environment for mergers is the disk of an
AGN (Bartos et al., 2017; Stone et al., 2017; McKernan et al., 2020; Tagawa et
al., 2020b). Theory and observations (Artymowicz et al., 1993; Levin &
Beloborodov, 2003; Tagawa et al., 2020b; Fan & Wu, 2022) suggest that black
holes get embedded in AGN disks due to capture via dynamical interactions
between the nuclear star cluster and the AGN disk (Ostriker, 1983; Miralda-
Escudé & Kollmeier, 2005) and by in-situ star formation (Levin & Beloborodov,
2003; Milosavljević & Loeb, 2004). The AGN disks then act as black hole
assembly lines (Cheng & Wang, 1999), bringing the black holes closer together
and helping them merge over relatively short time scales. It is also possible
that multiple black holes merge consecutively, resulting in unusually heavy
remnants (Yang et al., 2020; Tagawa et al., 2021a). Comparisons to the
observed black hole masses, spins and merger rate indicate that a sizable
fraction of the observed mergers may originate in AGNs (Gayathri et al.,
2021a; Tagawa et al., 2021a). AGNs could also explain some of the peculiar
detections, such as those with a high mass (Tagawa et al., 2021a) and possible
high eccentricity (Samsing et al., 2022; Tagawa et al., 2021b; Romero-Shaw et
al., 2020; Gayathri et al., 2022; Romero-Shaw et al., 2022).

Due to the gas-rich merger environment, a key signature of the AGN channel is
the possibility of electromagnetic emission accompanying the gravitational
wave 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 counterparts suggested so far.
These include the peculiar optical flare ZTF19abanrhr with luminosity $L\sim
10^{45}~{}{\rm{erg/s}}$ and duration $t_{\rm{duration}}\sim 1~{}\rm{month}$,
detected by the Zwicky Transient Facility (ZTF) in the AGN
${\rm{J124942.3+344929}}$. It was observed in spatial coincidence with, and
within $18~{}\rm{days}$ of GW190521 (Graham et al., 2020; Calderón Bustillo et
al., 2021), the heaviest black hole merger detected to date (Abbott et al.,
2020a). Recently six optical flares possibly associated with massive black
hole mergers with luminosity and duration similar to those of ZTF19abanrhr
have been additionally reported by Graham et al. (2022). The first
gravitational wave event, GW150914 (Abbott et al., 2016a), has also been
proposed to have an associated bright gamma-ray event detected by the Fermi
Gamma-ray Burst Monitor (GBM), GW150914-GBM (Connaughton et al., 2016, 2018)
with $L\sim 10^{49}~{}{\rm{erg/s}}$ and $t_{\rm{duration}}\sim 1~{}{\rm{s}}$
(but see Greiner et al. 2016; Savchenko et al. 2016). Finally, another gamma-
ray burst candidate has been associated with the tentative gravitational wave
event, LVT151012 (Bagoly et al., 2016), with properties similar to those of
GW150914-GBM. As these associations remain debated and controversial (Greiner
et al., 2016; Savchenko et al., 2016; Calderón Bustillo et al., 2021), it is
necessary to consider whether a self-consistent physical model can, in fact,
explain all of the properties of these tentative associations.

The possibility of electromagnetic emission from merging black holes in AGN
disks has been actively studied recently, mainly focusing on optical flares
(Graham et al., 2020; Kimura et al., 2021; Wang et al., 2021; McKernan et al.,
2019). Additionally, multiple scenarios have been proposed to produce emission
similar to GW150914-GBM in more generic galactic environments (e.g. Janiuk et
al., 2017; Perna et al., 2016). However, no model has been able to explain all
of the properties, including the luminosity, delay time, duration and color of
electromagnetic transients. Furthermore, ZTF19abanrhr began brightening only
approximately 18 days after the merger, which has not been physically
justified (see § 3.3.2 below).

Figure 1.—  Schematic pictures of the breakout emission. (a) Breakout emission
from the head of a jet launched from a solitary black hole embedded in an AGN
disk. From a solitary black hole, the breakout emission is produced
episodically, after replenishment of gas to the black hole (Tagawa et al.,
2022b). (b) The emission from shocks produced around a jet after a merger. The
jet direction is aligned with the black hole spin direction, which is
reoriented by the merger. As a result, the jet again collides with unshocked
AGN gas, producing breakout emission after the merger.

A promising process for the production of detectable electromagnetic emission
from AGN-assisted black hole mergers is the launch of a relativistic jet
through the Blandford-Znajek effect (Tagawa et al. 2022b; see Fig. 1 for a
schematic illustration). We evaluated the properties and observability of non-
thermal and thermal emission from shocks emerging around such jets. In the AGN
disk, a black hole is often surrounded by a circum-black hole disk (Tagawa et
al., 2022b). When the circum-black hole disk is advection dominated, as
realized in the slim disk model expected in this environment (Abramowicz et
al., 1988), a magnetically dominated state can be realized owing to the
accumulation of magnetic flux in the vicinity of the black hole (Cao 2011; the
magnetic field can be further magnified by an outflow; see e.g. Liska et al.
2020). In such cases, a jet is expected to be launched from a spinning black
hole via the Blandford-Znajek process (Blandford & Znajek, 1977).

As the jet collides with unshocked gas in the AGN disk, strong forward and
reverse shocks are formed in the disk and the jet head, respectively, and
shocked material surrounds the sides of the jet (cocoon). During the early
phases, photons in the shocked medium cannot escape from the system because
they are surrounded by an optically thick disk. As the shock approaches the
surface of the AGN disk, photons begin escaping from the system, leading to
luminous thermal and non-thermal emission (Fig. 1a). Note that solitary black
holes produce similar breakout emission, but it will be difficult to observe
because of lower brightness compared to that from merging black holes and the
low duty-cycle of the breakout emission (Tagawa et al., 2022a): the jet self-
regulates to an episodic behaviour by creating a cavity around the black hole,
with the breakout emission existing only for a small fraction of each cycle
(Kimura et al., 2021; Tagawa et al., 2022b). On the other hand, we show that
breakout emission from merging black holes is expected to be observable due to
its brightness, and that it accompanies up to several percent of black hole
mergers in these environments (see $\S$4.1). We further present a prescription
for how non-thermal and thermal emission from these shocks can be used to
”reverse engineer” the properties of the merger environments, applying the
technique to the transients proposed to be associated with black hole mergers:
ZTF19abanrhr, GW150914-GBM, and LVT151012-GBM.

##  2\. Model Description

In this section, we describe our model of accretion onto black holes, the
properties of thermal and non-thermal emission, and our numerical choices. The
various parameters of our model are summarized in Table 2 in Appendix A.
Readers not interested in the model details may skip directly to the next
section, describing our results.

###  2.1. Accretion onto black holes

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

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

where $\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., 2022), and $a_{\rm BH}$ is the
dimensionless spin of the black hole. Since the power of a shock emerged
around the jet and the luminosity of radiation emitted from the shock are
roughly proportional to the jet kinetic luminosity, the accretion rate onto
the black hole is a key quantity to determine the observed luminosity from the
system.

The accretion rate onto a circum-black hole disk in the AGN disk is evaluated
via the Bondi-Hoyle-Lyttleton rate, as given by Eq. (1) of Tagawa et al.
(2022b). By considering the reduction from the Bondi-Hoyle-Lyttleton rate, we
parameterized the fraction of the accretion rate onto the black hole
(${\dot{m}}$) over the Bondi-Hoyle-Lyttleton rate (${\dot{m}}_{\rm BHL}$) as
$f_{\rm acc}={\dot{m}}/{\dot{m}}_{\rm BHL}$. For example, low $f_{\rm acc}$
may be predicted due to winds from an accretion disk with a super-Eddington
rate, although recent simulations suggest that the conversion to wind is
moderate (Kitaki et al., 2021) for accretion flows onto black holes in an AGN
disk, in which the circularization radius is much larger than the trapping
radius. In addition, the accretion rate onto the black hole in a cavity in
active phases is estimated to be lower by a factor of a few compared to that
without a cavity (Tagawa et al., 2022b). On the other hand, for merger
remnants, inflow rates can be enhanced by shocks arising in the marginally
bound annuli of the post-merger circum-black hole disks, due to recoil kicks
imparted on the remnants. Considering this process, we adopt $f_{\rm acc}=15$
as a fiducial value, as discussed in detail in Appendix B.

Figure 2.— Various quantities as a function of the distance from the
supermassive black hole ($R$) for breakout emission. (a) The timescales
representing delay ($t_{\rm delay}$, solid black), breakout ($t_{\rm break}$,
dashed orange), and transparency ($t_{\rm trans}$, dotted blue). (b) The
observed duration of breakout emission ($t_{\rm duration}$, solid black), the
duration of the shock breakout ($t_{\rm breakout}$, dashed black), the angular
timescale ($t_{\rm ang}$, dotted orange), and the transparent timescale
($t_{\rm trans}$, dotted-dashed blue). (c) The shock luminosity (solid black),
the breakout luminosity by thermal emission (dashed black), and that by the
non-thermal emission ($L_{\rm non-thermal}$, orange). (d) The absorption
($\gamma_{\rm{a}}$, solid), minimum ($\gamma_{\rm{m}}$, dashed), and maximum
($\gamma_{\rm{max}}$, dotted) Lorentz factors. (e) The breakout temperature
(solid), the absorption temperature ($h\nu_{\rm{a}}/k_{\rm{B}}$, dashed), and
the minimum temperature ($h\nu_{\rm{m}}/k_{\rm{B}}$, dotted). (f) The shock
velocity $\gamma_{\rm FS}\beta_{\rm FS}$ (solid), $\gamma_{\rm FS}$ (dashed),
and $\beta_{\rm FS}$ (dotted). The black hole locations adopted in the
fiducial model are indicated with filled circles superposed on the black solid
lines. The dotted holizontal lines in panels (a)–(c) present the observed
values for $t_{\rm delay}$, $t_{\rm duration}$, and $L$ at $\sim
2~{}{\rm{eV}}$.

###  2.2. Thermal emission at breakout

Here we provide details on the mechanisms involved in the production of
thermal emission from shocks produced by the interaction between the AGN disk
and the jet launched from the accreting black holes. We consider emission from
the shocks propagating in optically thick media (an AGN disk, Fig. 1). We
assume that the black hole is at the midplane of the AGN disk, and the jet
direction is perpendicular to the AGN disk plane. On the other hand, if the
jet is inclined with respect to the angular momentum direction of the AGN disk
by an angle $i$, the distance to the surface of the AGN disk from the black
hole needs to be enhanced by a factor of $1/{\rm cos}(i)$ compared to the
estimates in this paper, while other modifications would be minor. While the
diffusion of photons is slower than the propagation of a shock ($v_{\rm
FS}=\beta_{\rm FS}c>c/\tau$), where $\tau$ is the optical depth of the AGN
disk, photons are trapped in the shock. Once the shock propagates to the
height at which the diffusion to the edge (the disk photosphere) becomes
faster than the propagation ($\beta_{\rm FS}c<c/\tau$), photons start to
escape. The observed properties (such as the temperature, duration, and
luminosity) of the emission from the shock breakout are differently
characterized depending on the ranges of the shock velocity (Fig. 2), as
described in the following (Ito et al., 2020a; Levinson & Nakar, 2020).

For Newtonian shocks with $\beta_{\rm FS}\gamma_{\rm FS}\lesssim 0.03$ (Nakar
& Sari, 2010; Sapir et al., 2013), where $\gamma_{\rm FS}$ is the Lorentz
factor of the forward shock, the radiation is in thermodynamic equilibrium as
there is time to emit a sufficient number of photons by free-free emission.
From the jump condition for a strong shock, the thermal equilibrium
temperature is

$\displaystyle T_{\rm breakout}$ $\displaystyle\sim(18\rho_{\rm AGN}v_{\rm
FS}^{2}/7a)^{1/4}$ $\displaystyle\sim 10^{4}\,{\rm K}\left(\frac{\beta_{\rm
FS}}{0.02}\right)^{1/2}\left(\frac{\rho_{\rm AGN}}{1\times 10^{-16}\,{\rm
g~{}cm^{-3}}}\right)^{1/4},$ (2)

where $a$ is the radiation constant. The breakout emission typically peaks in
the optical or UV bands, $E_{\rm\gamma,BO}\sim 2.8~{}k_{\rm{B}}T_{\rm
breakout}\sim 1$–$100$ eV, where $k_{\rm{B}}$ is the Boltzmann constant.

In non-relativistic regime, photons inside the shock start to diffuse out from
the AGN disk, when the photon diffusion time, $t_{\rm diff}\sim~{}d_{\rm
edge}^{2}\kappa_{\rm AGN}\rho_{\rm AGN}/c$, becomes equal to the shock
expansion time, $t_{\rm dyn}\sim~{}d_{\rm edge}/v_{\rm FS}$, where $d_{\rm
edge}$ is the thickness of the AGN disk above the shock, and $\kappa_{\rm
AGN}$ is the opacity of the AGN disk for thermal photons. By equating the
timescales, the thickness at the breakout is given by $d_{\rm
edge,BO}\sim~{}c/(v_{\rm FS}\kappa_{\rm AGN}\rho_{\rm AGN})$, and the duration
of emission from a breakout shell is

$\displaystyle t_{\rm breakout}$ $\displaystyle\sim c/(v_{\rm
FS}^{2}\kappa_{\rm AGN}\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},$ (3)

where we adopt the Thomson scattering opacity of $\kappa_{\rm AGN}\sim
0.4\,{\rm g/cm^{2}}$ assuming ionized gas. The observed duration for bright
emission ($t_{\rm duration}$) is typically given by $t_{\rm breakout}$ in non-
relativistic cases. The luminosity roughly evolves following
$L\propto{\rm{e}}^{t^{2}}$ (Sapir et al., 2011) and $L\propto t^{-4/3}$ before
and after breakout during a planar phase ($t<t_{\rm sph}$), respectively,
where $t_{\rm sph}=H_{\rm AGN}/v_{\rm FS}$ is the time of transition between
the planar and spherical geometries of the breakout shell after the breakout
of the shock. The dependence of the time evolution of the luminosity on the
density profile of the ambient material is found to be weak (Sapir et al.,
2011). After $t_{\rm sph}$, the luminosity evolves more modestly depending on
the pre-shocked density profile. Note that for simplicity we adopt formulae
obtained in a spherical geometry, while the evolution of the breakout
luminosity from shocks propagating in a disk geometry would require
modification from them (Grishin et al., 2021).

The delay time between the production of a jet and the breakout is roughly
given by

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

where the prefactor of $3/5$ accounts for the deceleration during the non-
relativistic regime. The time delay between the gravitational wave and the
breakout emission ($t_{\rm delay}$) is roughly given by $t_{\rm break}$ in
non-relativistic regimes. Note that $t_{\rm delay}$ can be lowered compared to
$t_{\rm break}$ by up to the breakout timescale of gravitational waves
($t_{\rm GW,break}\sim H_{\rm AGN}/c$) if the line-of-sight direction and the
jet propagating direction coincide. This is because the breakout emission is
produced closer to us in such cases compared to those for random jet
directions.

The luminosity of the shock breakout and the shock luminosity are typically
(Nakar & Sari, 2010)

$\displaystyle L_{\rm breakout}$ $\displaystyle\sim L_{\rm
sh}\sim\pi\theta_{\rm j}^{2}H_{\rm AGN}^{2}\rho_{\rm AGN}v_{\rm FS}^{3}$
$\displaystyle\sim 5\times 10^{43}\,{\rm erg/s}\left(\frac{\theta_{\rm
j}}{0.05}\right)^{2}$ $\displaystyle\left(\frac{H_{\rm AGN}}{5\times
10^{16}\,\mathrm{cm}}\right)^{2}\left(\frac{\rho_{\rm AGN}}{1\times
10^{-16}\,{\rm g~{}cm^{-3}}}\right)\left(\frac{\beta_{\rm
FS}}{0.1}\right)^{3},$ (5)

where $\theta_{\rm j}$ is the opening angle of the jet, which is given by
$\theta_{\rm j}\sim(L_{\rm j}\theta_{0}^{6}\beta_{\rm FS}^{2}/\rho_{\rm
AGN}H_{\rm AGN}^{2}c^{3})^{1/10}$ for non-relativistic regimes (Bromberg et
al., 2011).

For Newtonian shocks with $0.03\lesssim\beta_{\rm FS}\gamma_{\rm FS}\lesssim
1$, the breakout temperature ($T_{\rm breakout}$) strongly depends on the
shock velocity (Budnik et al., 2010; Sapir et al., 2013). In these regimes,
since the number of photons produced by free-free emission during propagation
is lower than that required to establish thermal equilibrium (Nakar & Sari,
2010), the radiation is out of thermodynamic equilibrium and photons roughly
follow local Compton equilibrium, whose radiation is characterized by a Wien
spectrum $de_{\rm ph}/d\nu\propto\nu^{3}{\rm{e}}^{-h\nu/k_{\rm{B}}T_{\rm
breakout}}$, where $\nu$ and $e_{\rm ph}$ are the frequency and the spectral
energy density of photons, and $h$ is the Plank constant. The breakout
temperature of the photons is calculated and fitted as (Sapir et al., 2013)

$\displaystyle{\rm log}_{10}\left(\frac{T_{\rm breakout}}{{\rm 10^{4}~{}{\rm
K}}}\right)=0.975+1.735\left(\frac{\beta_{\rm FS}}{0.1}\right)^{1/2}$
$\displaystyle+\left[0.26-0.08\left(\frac{\beta_{\rm
FS}}{0.1}\right)^{1/2}\right]{\rm log}_{10}\left(\frac{n_{\rm
AGN}}{10^{15}~{}{\rm cm^{-3}}}\right),$ (6)

where $n_{\rm AGN}=\rho_{\rm AGN}/m_{\rm{p}}$ is the number density of the AGN
disk and $m_{\rm{p}}$ is the proton mass. Note that this formula likely
somewhat underestimates $T_{\rm breakout}$ for high $\beta_{\rm FS}$ and flat
density profiles (roughly up to a factor of a few, Fig. 3 in Sapir et al.
(2013)).

The duration of the breakout emission ($t_{\rm duration}$), the breakout
timescale ($t_{\rm break}$), and the shock luminosity ($L_{\rm sh}$) in the
fast-Newtonian regimes are given by Eqs. (2.2), (2.2), and (2.2) as in the
slow-Newtonian regimes.

We follow Nakar & Sari (2012) to model evolution in relativistic regimes. For
relativistic shocks with $\beta_{\rm FS}\gamma_{\rm FS}\gtrsim 1$, pairs are
abundantly created, which enhances the production of photons inside the shock
and regulates the temperature in the rest frame of the downstream plasma to
$\sim 100-200\,{\rm keV}$, almost independent of $\beta_{\rm FS}$ (Budnik et
al., 2010). Here, the relativistic shells are accelerated from the initial
Lorentz factor of $\gamma_{\rm FS}$ to $\gamma_{\rm FS,f}=\gamma_{\rm
FS}^{1+\sqrt{3}}$ due to the decrease of the gas density before breakout for
$\gamma_{\rm FS}\lesssim 4$ (Johnson & McKee, 1971), which is satisfied below.
Note that breakout emission for $\gamma_{\rm FS}\gtrsim 4$ is not well
investigated and applicability of the formalism is unclear. We also assume
that the observer is within the angle of $\lesssim 1/\gamma_{\rm sf,f}$ from
the direction of the shock propagation, and that radiation is beamed, unless
stated otherwise, where $\gamma_{\rm sf,f}\approx\gamma_{\rm sf}^{1+\sqrt{3}}$
is the final Lorentz factor of the shocked fluid, $\gamma_{\rm sf}$ is the
initial Lorentz factor of the shocked fluid, and $\gamma_{\rm
FS}\sim\sqrt{2}\gamma_{\rm sf}$ for strong shocks in relativistic regimes.

The delay time between a jet launch and the breakout from the AGN disk is
roughly given by

$\displaystyle t_{\rm delay}\sim{\rm{max}}\left(t_{\rm break,rel},t_{\rm
trans}\right),$ (7)

where

$\displaystyle t_{\rm break,rel}=\frac{H_{\rm AGN}}{4\gamma_{\rm FS}^{2}c}$
(8)

is the breakout timescale of the shock in relativistic regimes, in which the
factor of $1/2$ accounts for the deceleration during Blandford-McKee phases
(Blandford & McKee, 1976) and the factor of $H_{\rm AGN}/(2\gamma_{\rm
FS}^{2}c)$ approximately accounts for the difference between the apparent
travel time of the gravitational waves and the forward shock ($\sim~{}H_{\rm
AGN}/(\beta_{\rm FS}c)-H_{\rm AGN}/c$). $t_{\rm trans}$ is the time at which
the breakout shell with the temperature being $k_{\rm{B}}T_{\rm breakout}\sim
200~{}{\rm keV}~{}\gamma_{\rm sf,f}$ in the shock immediate downstream becomes
optically thin by the annihilation of pairs at the temperature of

$\displaystyle k_{\rm{B}}T_{\rm breakout}\sim 50~{}{\rm keV}~{}\gamma_{\rm
sf,f}$ (9)

due to cooling of the shell by adiabatic expansion (we adopt Eq. 13 of Nakar &
Sari (2012)). The timescale $t_{\rm trans}$ should be taken into account for
relativistic regimes since shocks are not transparent before the annihilation
of pairs due to the enhancement of the optical depth by a factor of $\gtrsim
100$.

The duration of the emission from the breakout shell is

$\displaystyle t_{\rm duration}\sim{\rm{max}}\left(t_{\rm ang},t_{\rm
trans},t_{\rm breakout}\right),$ (10)

where

$\displaystyle t_{\rm ang}=\frac{H_{\rm AGN}}{2\gamma_{\rm sf,f}^{2}c}$ (11)

is the timescale during which radiation emitted by the same shell is observed
due to its round shape in the relativistic regime. In relativistic regimes and
when shocks are propagating towards the observer, the breakout timescale
($t_{\rm breakout}$) is reduced due to the beaming effect by $\sim\gamma_{\rm
FS}^{2}$, pair production by $100~{}\gamma_{\rm FS}$, and Klein-Nishina effect
by $\gamma_{\rm FS}^{-2}$ (Budnik et al., 2010) as

$\displaystyle t_{\rm breakout}\sim\frac{c}{100v_{\rm FS}^{2}\kappa_{\rm
AGN}\rho_{\rm AGN}\gamma_{\rm FS}}.$ (12)

For relativistic cases, the shock luminosity is (Blandford & McKee, 1976)

$\displaystyle L_{\rm sh}\sim\pi\theta_{\rm j}^{2}H_{\rm AGN}^{2}\rho_{\rm
AGN}v_{\rm FS}\gamma_{\rm FS}^{2}c^{2},$ (13)

while the isotropic equivalent breakout luminosity is given by (Nakar & Sari,
2012)

$\displaystyle L_{\rm breakout}\sim\pi f_{\rm beaming}\theta_{\rm j}^{2}H_{\rm
AGN}^{2}\rho_{\rm AGN}v_{\rm FS}\gamma_{\rm FS}^{2}c^{2}\frac{\gamma_{\rm
sf,f}}{4\gamma_{\rm sf}}$ (14)

where $f_{\rm beaming}$ is the beaming factor, takes into account the angular
effect due to the fact that emission is concentrated in the direction of shock
motion. We assume $f_{\rm beaming}=2\gamma_{\rm sf,f}^{2}$ approximating that
the radiation is concentrated within the angle $\sim 1/\gamma_{\rm sf,f}$ from
the direction of the shock motion. Eq. (13) includes an additional boosting
factor to account for conversion of the emitted power to the observer frame,
to obtain the full Doppler boost. We assume that the breakout shell has a
width of $v_{\rm FS}t_{\rm duration}$ (which is different from the width for
the non-thermal emission of $v_{\rm FS}t_{\rm breakout}$). The factor of
$\gamma_{\rm sf,f}/4\gamma_{\rm sf}$ in Eq. (14) represents the reduction of
the internal energy due to adiabatic expansion, since the breakout shell
becomes transparent at the temperature of $\sim 50~{}{\rm keV}~{}\gamma_{\rm
sf,f}$, which is lower than the temperature of the shell in the immediately
shock downstream ($\sim 200~{}{\rm keV}~{}\gamma_{\rm sf}$) by the factor.

In summary, we adopt Eqs. (2.2), (2.2), (9) for $T_{\rm breakout}$ for
$\beta_{\rm FS}\gamma_{\rm FS}\leq 0.03$, $1\geq\beta_{\rm FS}\gamma_{\rm
FS}\geq 0.03$, and $\beta_{\rm FS}\gamma_{\rm FS}\geq 1$, respectively. The
radiation is roughly characterized by a black-body radiation and a Wien
spectrum for $\beta_{\rm FS}\gamma_{\rm FS}\leq 0.03$ and $\beta_{\rm
FS}\gamma_{\rm FS}\geq 0.03$, respectively, although numerical simulations
find that the spectrum at breakout is typically softer than that, especially
for Newtonian slow shocks and in relativistic regimes (Ito et al., 2020a;
Budnik et al., 2010). Also, we assume that $t_{\rm delay}$ is given by Eqs.
(2.2) and (7), $t_{\rm duration}$ is given by Eqs. (2.2) and (10) for
$\beta_{\rm FS}\gamma_{\rm FS}\leq 1$ and $\beta_{\rm FS}\gamma_{\rm FS}\geq
1$, respectively, and $L_{\rm sh}$ and $L_{\rm breakout}$ are given by Eqs.
(2.2) for $\beta_{\rm FS}\gamma_{\rm FS}\leq 1$, and are, respectively, given
by Eqs. (13) and (14) for $\beta_{\rm FS}\gamma_{\rm FS}\geq 1$.

The velocity of the shocked fluid ($\beta_{\rm sf}$) is determined by the
competition between the power of the jet $L_{\rm j}$ and the mass of the
ambient material swept up by the jet (Bromberg et al., 2011), which is derived
using the density ($\rho_{\rm AGN}$), the scale height of the AGN disk
($H_{\rm AGN}$), and the injected opening angle of the jet ($\theta_{\rm j}$).
Note that at the transitions where $\beta_{\rm sf}\gamma_{\rm sf}\sim 1$ or
${\tilde{L}}\sim 1$, where ${\tilde{L}}$ is the ratio between the energy
density of the jet and the rest-mass energy density of the surrounding medium
at the location of the head, the formulae in Bromberg et al. (2011)
unexpectedly predict that $\beta_{\rm sf}\gamma_{\rm sf}$ decreases as $L_{\rm
j}$ increases. To avoid this, we set the upper limit for ${\tilde{L}}\leq 1$
and the lower limit for ${\tilde{L}}\geq 1$ to be $\beta_{\rm sf}\gamma_{\rm
sf}=1$, since the two cases correspond to $\beta_{\rm sf}<1$ (with
$\gamma_{\rm sf}\sim 1$) and $\gamma_{\rm sf}>1$ (with $\beta_{\rm sf}\sim
1$), respectively. Thus, $\beta_{\rm sf}\gamma_{\rm sf}$ is not accurately
determined at around $\beta_{\rm sf}\gamma_{\rm sf}\sim 1$.

To estimate the properties of breakout emission from the shock, we assume that
the breakout velocity is given by the head velocity of the shocked region,
although the directions of the jets may be random (Tagawa et al., 2020a) and
the side of the shock (cocoon) may first break out. Since the head velocity is
$\sim 5(\theta_{\rm c}/0.2)^{-1}$ times faster than the cocoon expansion
velocity (Bromberg et al., 2011), where $\theta_{\rm c}$ is the opening angle
of the cocoon, we assume that the cocoon can break out faster than the jet
head when the angle between the propagation direction of the jet head and the
AGN disk plane is less than $\theta_{\rm c}$. Assuming that the jet direction
is random, the probability that the shock breaks out the AGN surface from its
head is $\sim\int_{0}^{\pi/2-\theta_{\rm c}}{\rm
sin}(\theta)d\theta/\int_{0}^{\pi/2}{\rm sin}(\theta)d\theta\sim 1-\theta_{\rm
c}\sim 0.8$ for $\theta_{\rm c}=0.2$. Thus, the shock emission most likely
breaks out the AGN disk surface from the head, and we only consider this case.

###  2.3. Non-thermal emission

In this subsection, we summarize the properties of the non-thermal emission
emerged after the shock breakout. While photons are trapped in the shock
($\beta_{\rm FS}c>c/\tau$), the shock is mediated by scattering with the
photons, in which electrons are not (or inefficiently) accelerated (Ito et
al., 2020a). Once photons start to escape, a collisionless shock can be
formed, where electrons can be accelerated (e.g. Kashiyama et al., 2013; Ito
et al., 2020b). These electrons produce non-thermal emissions via synchrotron
and inverse Compton scattering processes. Here we focus on non-thermal
emission after the shock breakout timescale, which is given by $t_{\rm break}$
(Eq. 2.2) for non-relativistic and $t_{\rm break,rel}$ (Eq. 7) for
relativistic regimes, respectively. We focus on the non-thermal emission
around the peak of the lightcurve. The quantitative estimate of the time
evolution would require numerical simulations, which is beyond the scope of
this paper.

After the breakout of the shock, the non-thermal luminosity decreases
depending on the density profile of the AGN disk and the speed of the shock.
If the non-thermal radiation emerges – and hence it is observed – over a
timescale longer than the intrinsic duration of particle acceleration, its
brightness will be correspondingly diluted by a factor which we indicate with
$f_{\rm dilution}$. For non-relativistic shocks ($\beta_{\rm FS}\gamma_{\rm
FS}\leq 1$), we have $f_{\rm dilution}=1$ since the two timescales are
identical. For relativistic shocks in gas-pressure dominated regimes for the
AGN disk, since the gas density gradually changes as a function of the
distance from the AGN disk plane, the luminosity gradually decreases on the
observational timescale, and hence we again assume $f_{\rm dilution}=1$. On
the other hand, for relativistic shocks in radiation-pressure dominated
regimes for the AGN disk (which is realized for $R\lesssim 2\times
10^{-2}~{}{\rm pc}$ in the fiducial model, e.g. Haiman et al. 2009), the gas
density abruptly decreases at a few $H_{\rm AGN}$ (e.g. Grishin et al., 2021).
In such cases, the non-thermal luminosity also abruptly decreases and is high
only within the timescale of $\sim t_{\rm breakout}$. In this case the
observed emission is diluted on the timescale $t_{\rm dur,NT}={\rm max}(t_{\rm
ang},t_{\rm breakout})$. By assuming for simplicity that the gas density above
the AGN disk is almost in a vacuum, we then set $f_{\rm dilution}=t_{\rm
breakout}/t_{\rm dur,NT}$ for the cases with relativistic shocks in radiation-
pressure dominated regimes. We should note that the assumption may
underestimate the non-thermal emission, because disk winds might exist that
increase the value of $f_{\rm dilution}$.

We assume that a fraction $\epsilon_{\rm{e}}$ of the kinetic energy of the
shock is used to accelerate electrons in a collisionless shock (Sironi et al.,
2013). Then, the synchrotron luminosity is given by (Fan & Piran, 2008)

$\displaystyle L_{\rm syn}$ $\displaystyle\sim$ $\displaystyle\frac{L_{\rm
non-thermal}}{1+Y_{\rm SSC}+Y_{\rm 2ndIC}}$ (15) $\displaystyle\sim$
$\displaystyle\epsilon_{\rm{e}}f_{\rm beaming}f_{\rm dilution}\frac{L_{\rm
sh}}{1+Y_{\rm SSC}+Y_{\rm 2ndIC}}$

where $Y_{\rm SSC}$ and $Y_{\rm 2ndIC}$ are the powers of synchrotron-self
Compton and second order inverse Compton scattering compared to that of
synchrotron emission, and $L_{\rm non-thermal}$ is the luminosity by non-
thermal emission. We ignore the inverse Compton scattering of the thermal
photons of the AGN disk and those of breakout emission. The former is because
its energy density is much lower compared to that of the synchrotron photons.
The latter is because the scattering of such photons are typically in the
Klein-Nishina regimes in our case. For $\gamma_{\rm sf}\beta_{\rm sf}\leq 1$,
we set $f_{\rm beaming}=1$, and otherwise $f_{\rm beaming}=2\gamma_{\rm
sf,f}^{2}$.

Here, plasma and/or MHD instabilities are assumed to amplify the magnetic
field to $\epsilon_{\rm{B}}\lesssim 10^{-3}$–$10^{-1}$, while electrons are
accelerated via the first-order Fermi process with an energy fraction of
$\epsilon_{\rm{e}}\lesssim 10^{-2}$–$0.3$ (from observations, e.g., Panaitescu
& Kumar 2001; Frail et al. 2005, and from theoretical studies, e.g.,
Spitkovsky 2008; Sironi et al. 2013). We assume that electrons are accelerated
in the shock to a power-law distribution of Lorentz factors
$\gamma^{\prime}_{\rm{e}}$ as
$N(\gamma_{\rm{e}}^{\prime})d\gamma_{\rm{e}}^{\prime}\propto{\gamma_{\rm{e}}^{\prime}}^{-p}d\gamma_{\rm{e}}^{\prime}$
with a minimum ($\gamma_{\rm{m}}^{\prime}$) and a maximum value
($\gamma_{\rm{max}}^{\prime}$) and a power-law index $p$, where primes are
used for quantities in the fluid comoving frame. Assuming that a fraction
$\epsilon_{\rm{e}}$ of the kinetic energy is converted to electron energy, the
minimum Lorentz factor $\gamma_{\rm{m}}^{\prime}$ is given by

$\displaystyle\gamma_{\rm{m}}^{\prime}\sim\epsilon_{\rm{e}}\left(\frac{p-2}{p-1}\right)\frac{m_{\rm{p}}}{m_{\rm{e}}}(\gamma_{\rm
sf}-1)\sim
40~{}\left(\frac{\epsilon_{\rm{e}}}{0.3}\right)\left(\frac{\gamma_{\rm
sf}-1}{0.2}\right),$ (16)

where $m_{\rm{e}}$ is the electron mass.

By comparing the cooling timescale by the synchrotron radiation and the
acceleration timescale by the first-order Fermi acceleration mechanism, the
maximum Lorentz factor of electrons is determined as

$\displaystyle\gamma_{\rm{max}}^{\prime}$ $\displaystyle=\frac{v_{\rm
sf}}{c}\left(\frac{9\pi e}{10\sigma_{\rm T}B^{\prime}_{\rm
sf}\xi}\right)^{1/2}$ $\displaystyle\sim 1\times
10^{6}\xi^{-1/2}\left(\frac{\beta_{\rm
sf}}{0.5}\right)\left(\frac{A_{\gamma_{\rm sf}}}{1.5}\right)^{-1/4}$
$\displaystyle\left(\frac{\epsilon_{\rm{B}}}{0.1}\right)^{-1/4}\left(\frac{\rho_{\rm
AGN}}{1\times 10^{-16}\,{\rm g~{}cm^{-3}}}\right)^{-1/4}$ (17)

where $e$ is the electric charge, $A_{\gamma_{\rm sf}}=(\gamma_{\rm
sf}-1)(4\gamma_{\rm sf}+3)$,

$\displaystyle B^{\prime}_{\rm sf}$
$\displaystyle=(8\pi\epsilon_{\rm{B}}e^{\prime}_{\rm sf})^{1/2}$
$\displaystyle\sim 6\times 10^{2}~{}{\rm G}~{}\left(\frac{A_{\gamma_{\rm
sf}}}{1.5}\right)^{1/2}\left(\frac{\epsilon_{\rm{B}}}{0.1}\right)^{1/2}\left(\frac{\rho_{\rm
AGN}}{1\times 10^{-16}\,{\rm g~{}cm^{-3}}}\right)^{1/2},$ (18)

is the magnetic field in the shocked medium (parameterized in terms of the
fraction $\epsilon_{\rm{B}}$ of postshock energy carried by the post shock
magnetic field), $e^{\prime}_{\rm sf}=A_{\gamma_{\rm sf}}\rho_{\rm AGN}c^{2}$
is the internal energy density of the shocked medium, and $\xi$ is the
parameter representing the ratio of the mean free path to the Larmor radius of
electrons. We use $\xi=1$, but this choice only affects $\nu_{\rm{max}}$ and
has little impact on most of the conclusions in this paper.

The cooling Lorentz factor of electrons, with which electrons can cool in the
dynamical timescale due to radiation, is given by

$\displaystyle\gamma^{\prime}_{\rm c}={\rm{max}}\left(1,\frac{6\pi
m_{\rm{e}}c}{\sigma_{\rm T}B_{\rm sf}^{\prime 2}\gamma_{\rm sf}t_{\rm
break}}\right)$ (19)

assuming that inverse Compton scattering is subdominant for cooling. Note that
non-thermal emission is characterized by fast cooling regimes
($\gamma^{\prime}_{\rm{m}}>\gamma^{\prime}_{\rm c}$) in the fiducial model.
The cooling timescale for $\gamma_{\rm{m}}^{\prime}$ is

$\displaystyle t_{\rm c}(\gamma_{\rm{m}}^{\prime})\sim 3\times
10\,{\rm{s}}~{}\left(\frac{\epsilon_{\rm{B}}}{0.1}\right)^{-1}\left[\frac{(\gamma_{\rm
sf}-1)(4\gamma_{\rm sf}+3)}{1.5}\right]^{-1}$
$\displaystyle\left(\frac{\rho_{\rm AGN}}{1\times 10^{-16}\,{\rm
g~{}cm^{-3}}}\right)^{-1}\left(\frac{\gamma_{\rm{m}}^{\prime}}{40}\right)^{-1}\left(\frac{\gamma_{\rm
sf}^{\prime}}{1.2}\right)^{-1}.$ (20)

From $t_{\rm c}(\gamma_{\rm{m}}^{\prime})$, the typical shell width of
electrons with $\gamma_{\rm{m}}^{\prime}$ emitting the synchrotron photons is
approximated as $\Delta_{\rm shell}(\gamma_{\rm{m}}^{\prime})\sim t_{\rm
c}(\gamma_{\rm{m}}^{\prime})v_{\rm FS}\sim 5\times 10^{11}\,{\rm cm}~{}[t_{\rm
c}(\gamma_{\rm{m}}^{\prime})/30~{}{\rm{s}}](\beta_{\rm FS}/0.5)$.

The Lorentz factor below which synchrotron self-absorption is effective is
(Rybicki & Lightman, 1979; Fouka & Ouichaoui, 2011)

$\displaystyle\gamma_{\rm{a}}^{\prime}=\gamma_{\rm{m}}^{\prime}\times\left(\tau_{\rm
q}C_{q+1}\right)^{1/(q+4)},$ (21)

where

$\displaystyle\tau_{q}=\frac{\pi}{3\sqrt{3}}\frac{(q^{2}+q-2)\gamma_{\rm
1}^{\prime-5}}{1-(\gamma_{\rm 2}^{\prime}/\gamma_{\rm
1}^{\prime})^{-q+1}}\frac{e\rho^{\prime}_{\rm sf}\Delta_{\rm
shell}(\gamma_{\rm{a}}^{\prime})}{B^{\prime}_{\rm sf}{\rm sin}\theta_{\rm
PA}m_{\rm{p}}},$ (22)

$\theta_{\rm PA}$ is the typical pitch angle between the magnetic field and
the velocity of electrons, $\gamma_{2}$ and $\gamma_{1}$ are the maximum and
minimum Lorentz factors of non-thermal electrons with power-law index of $q$
around $\gamma_{\rm{a}}$ while the self-absorption heats electrons,
respectively,

$\displaystyle
C_{q}=\frac{2^{(q+1)/2}}{q+1}\Gamma\left(\frac{q}{4}-\frac{1}{12}\right)\Gamma\left(\frac{q}{4}+\frac{19}{12}\right),$
(23)

$\Gamma$ is the Gamma function, $C_{p+1}\sim 1.6$ for $q=2.5$,
$n^{\prime}_{\rm norm}$ is the normalization for the electron number density,
and we assume that electrons are randomly oriented in the frame of the shocked
fluid. We set $\gamma_{2}^{\prime}=\gamma_{\rm{max}}^{\prime}$,
$\gamma_{1}^{\prime}=\gamma_{\rm{m}}^{\prime}$, and $q=p$ for
$\gamma_{\rm{m}}^{\prime}<\gamma_{\rm{a}}^{\prime}$, and
$\gamma_{2}^{\prime}=\gamma^{\prime}_{\rm{m}}$,
$\gamma_{1}^{\prime}=\gamma_{\rm{a}}^{\prime}$, and $q=2$ for
$\gamma_{\rm{a}}^{\prime}<\gamma_{\rm{m}}^{\prime}$. We determine
$\gamma_{\rm{a}}^{\prime}$ to satisfy Eqs. (21) and (22). For
$\gamma_{\rm{m}}^{\prime}<\gamma_{\rm{a}}^{\prime}$, we assume that the power-
law slope of the Lorentz factor of electrons during synchrotron self-
absorption is the same as the one injected ($p$) at around
$\gamma^{\prime}=\gamma_{\rm{a}}^{\prime}$. If synchrotron self-absorption
becomes effective after the synchrotron cooling modifies the electron
distribution to the slope of $p-1$, $\gamma_{\rm{a}}^{\prime}$ is enhanced by
a factor of $\sim 1.08$ in the fiducial model (Eq. 21).

With the synchrotron frequency given by

$\displaystyle\nu_{\rm sync}$ $\displaystyle=\gamma_{\rm
sf}\frac{eB^{\prime}_{\rm sf}}{2\pi m_{\rm{e}}c}$ $\displaystyle\simeq 2\times
10^{9}~{}{\rm Hz}\left(\frac{\gamma_{\rm
sf}}{1.2}\right)\left(\frac{B^{\prime}_{\rm sf}}{600~{}{\rm G}}\right),$ (24)

the frequencies corresponding to $\gamma_{\rm{m}}^{\prime}$,
$\gamma_{\rm{max}}^{\prime}$, and $\gamma_{\rm{a}}^{\prime}$ are,
respectively, (Fan & Piran, 2008)

$\displaystyle\nu_{\rm{m}}\simeq 3\times 10^{12}~{}{\rm
Hz}\left(\frac{\gamma_{\rm{m}}^{\prime}}{40}\right)^{2}\left(\frac{\nu_{\rm
sync}}{2\times 10^{9}~{}{\rm Hz}}\right),$ (25) $\displaystyle\nu_{\rm{max}}$
$\displaystyle\simeq 2\times 10^{21}~{}{\rm
Hz}\left(\frac{\gamma_{\rm{max}}^{\prime}}{1\times
10^{6}}\right)^{2}\left(\frac{\nu_{\rm sync}}{2\times 10^{9}~{}{\rm
Hz}}\right)$ $\displaystyle\simeq 2\times 10^{21}~{}{\rm
Hz}\left(\frac{\xi}{1}\right)^{-1}\left(\frac{\beta_{\rm sf}^{2}\gamma_{\rm
sf}}{0.3}\right),$ (26)

and

$\displaystyle\nu_{\rm{a}}=7\times 10^{12}~{}{\rm
Hz}\left(\frac{\gamma_{\rm{a}}^{\prime}}{60}\right)^{2}\left(\frac{\nu_{\rm
sync}}{2\times 10^{9}~{}{\rm Hz}}\right).$ (27)

For $\nu_{\rm{a}}<\nu_{\rm{m}}$, the luminosity by synchrotron radiation is
approximately given by (Fan & Piran, 2008)

$\displaystyle\nu L(\nu)\sim L_{\rm
syn}\times\left\\{\begin{array}[]{cl}(\nu/\nu_{\rm{m}})^{\frac{-p+2}{2}}~{}~{}~{}\mathrm{for}~{}~{}\nu_{\rm{m}}<\nu<\nu_{\rm{max}}\\\
(\nu/\nu_{\rm{m}})^{1/2}~{}~{}~{}\mathrm{for}~{}~{}\nu_{\rm{a}}<\nu<\nu_{\rm{m}}\\\
(\nu/\nu_{\rm{a}})^{3}(\nu_{\rm{a}}/\nu_{\rm{m}})^{1/2}~{}~{}~{}\mathrm{for}~{}~{}\nu_{\rm
syn}<\nu<\nu_{\rm{a}}.\\\ \end{array}\right.$ (31)

On the other hand, for $\nu_{\rm{m}}<\nu_{\rm{a}}$, electrons are thermalized
below and around $\gamma_{\rm{a}}$. We assume that the emission by synchrotron
radiation in such cases is approximately given by

$\displaystyle\nu L(\nu)\sim L_{\rm
syn}\times\left\\{\begin{array}[]{cl}(\nu/\nu_{\rm{m}})^{\frac{-p+2}{2}}+{\rm
exp}(1-(\nu/\nu_{\rm{a}})^{1/2})\\\
~{}~{}~{}\mathrm{for}~{}~{}\nu_{\rm{a}}<\nu<\nu_{\rm{max}}\\\
(\nu/\nu_{\rm{a}})^{3}~{}~{}~{}\mathrm{for}~{}~{}\nu<\nu_{\rm{a}},\end{array}\right.$
(34)

where the second term for $\nu_{\rm{a}}<\nu<\nu_{\rm{max}}$ and the luminosity
for $\nu<\nu_{\rm{a}}$ represent emission from thermalized electrons. The
order of $\nu_{\rm c}<\nu_{\rm{m}}<\nu_{\rm{a}}$ is astrophysically rare (Gao
et al., 2013), while it is often realized in the systems we investigated. This
is presumably because the ambient material is dense, and the outflow is
characterized by its high energy and compact size. The synchrotron-self
Compton spectrum of synchrotron photons by thermalized electrons is assumed to
be

$\displaystyle\nu L(\nu)\sim L_{\rm syn}Y_{\rm SSC}$
$\displaystyle\times\left\\{\begin{array}[]{cl}{\rm exp}(1-(\nu/\nu_{\rm
a,SSC})^{1/2})+(\nu/\nu_{\rm m,SSC})^{(2-p)/2}\\\
~{}~{}~{}\mathrm{for}~{}~{}\nu_{\rm a,SSC}<\nu\\\ (\nu/\nu_{\rm a,SSC})^{2}\\\
~{}~{}~{}\mathrm{for}~{}~{}\nu<\nu_{\rm a,SSC},\end{array}\right.$ (37)

whose thermal components are derived by considering inverse Compton scattering
of synchrotron photons by electrons thermalized by synchrotron self-absorption
(with the distribution $N_{\rm{e}}\propto\gamma^{2}{\rm
exp}(-\gamma/\gamma_{\rm{a}})$) (Ghisellini et al., 1988; Gao et al., 2013),
where $\nu_{\rm a,SSC}\approx(4/3)\gamma_{\rm{a}}^{\prime 2}\nu_{\rm{a}}$ and
$\nu_{\rm m,SSC}\approx(4/3)\gamma_{\rm{m}}^{\prime 2}\nu_{\rm{m}}$ are the
absorption and minimum frequencies for synchrotron-self Compton. Note that
since the low-energy range of the synchrotron-self Compton spectrum is
determined by the up-scattered photon distribution, the synchrotron-self
Compton spectrum is softer than the synchrotron spectrum. The spectrum of 2nd
order inverse Compton scattering is computed in a similar fashion as that of
the synchrotron-self Compton. For $\epsilon_{\rm{B}}\ll\epsilon_{\rm{e}}$,
$Y_{\rm SSC}$ and $Y_{\rm 2ndIC}$ are approximated as $Y_{\rm
SSC}=(\epsilon_{\rm{e}}/\epsilon_{\rm{B}})^{1/3}$ and $Y_{\rm
2ndIC}=(\epsilon_{\rm{e}}/\epsilon_{\rm{B}})^{2/3}$. Note that second order
inverse Compton scattering is in the Thomson regime if
${\rm{max}}(\nu_{\rm{m}}\gamma^{\prime 3}_{\rm{m}},\nu_{\rm{a}}\gamma^{\prime
3}_{\rm{a}})<m_{\rm{e}}c^{2}/2$ for the fast cooling regime (Fan & Piran,
2008). We do not consider the third order inverse Compton scattering since it
is suppressed by the Klein Nishina effect, where in the Klein-Nishina regimes
the scattering cross section is reduced compared to Thomson scattering due to
quantum electrodynamical corrections (Blumenthal & Gould, 1970).

Fig. 2 shows the dependence of the properties of emission on the distance $R$
from the supermassive black hole. It can be seen that this dependence is very
strong, which can be explained as follows. At smaller radii, the scale height
of the disk is lower and the density of the AGN disk is higher. Due to the
higher gas density, the accretion rate onto stellar-mass BHs is higher,
resulting in a stronger jet power and a higher shock velocity. The delay
timescale is also shorter for smaller scale height and higher shock velocity
(Eq. 2.2). Also, the duration is shorter for a higher AGN density and higher
shock velocity (Eq. 2.2). The temperature of the thermal emission is higher
for higher shock velocity and AGN density (Eq. 2.2). These are the reasons
that properties of emission strongly depend on the radius.

###  2.4. Numerical choices

When we discuss the association of ZTF19abanrhr to GW190521 and predict future
electromagnetic counterparts in Results, we use the following values for our
model parameters. As constrained by Graham et al. (2020) and Abbott et al.
(2020a), we assume that the masses of the merged black hole and the
supermassive black hole hosting the AGN are $m=150\,{\rm M}_{\odot}$ and
$M=10^{8}\,{\rm M}_{\odot}$, respectively. The jet energy conversion
efficiency to $\eta_{\rm j}=0.5$ (Tchekhovskoy et al., 2010), reflecting the
spin magnitude of the merged black hole ($a_{\rm BH,rem}\sim 0.7$ (Abbott et
al., 2020a)). In order to reproduce the observed luminosity of the AGN in
units of the Eddington Luminosity ($L_{\rm Edd}$) to be $\sim 0.2$ (Graham et
al., 2020) with a radiation efficiency of $\eta_{\rm rad}=0.1$, we set the gas
inflow rate from the outer boundary ($R_{\rm out}=10\,{\rm{pc}}$) of the AGN
disk to ${\dot{M}}_{\rm in}=50~{}{L}_{\rm Edd}/c^{2}$, the angular momentum
transfer parameter of the AGN disk to $m_{\rm AM}=0.6$ (Thompson et al.,
2005), and the viscous parameter of the AGN disk to $\alpha_{\rm AGN}=0.1$
(Martin et al., 2019). We assume that the opening angle of the injected jet is
$\theta_{0}=0.3$ (Hada et al., 2018; Berger, 2014), and the fraction of the
accretion rate onto the black hole (${\dot{m}}$) over the capture rate to
$f_{\rm acc}=15$ as discussed above. We adopt optimistic values for the
fraction of postshock energy carried by the post shock magnetic field and that
by electrons to $\epsilon_{\rm{B}}=0.1$ and $\epsilon_{\rm{e}}=0.3$ (e.g.,
Panaitescu & Kumar 2001; Frail et al. 2005; Sironi et al. 2013), respectively,
and the power-law slope for injected electrons accelerated by the first order
Fermi process to $p=2.5$.

##  3\. Results

###  3.1. Breakout emission from merging black holes

Breakout emission is expected in association with black hole mergers as
follows. The jet direction is aligned to the black hole spin direction. The
black hole spin directions before the mergers tend to be aligned perpendicular
to the AGN disk plane. Mergers are generally expected to occur with a
different spin direction since the angular momentum directions of merging
binaries are predicted to be quasi-randomized (Samsing et al., 2022) due to
frequent binary-single interactions (Tagawa et al., 2020a) in addition to
inhomogeneities in the outer regions of the AGN disk (Tagawa et al., 2021b).
This results in the reorientation of the jets. After this reorientation, as
the jet once again collides with unshocked gas, shocks emerge and breakout
emission is released following the emission of gravitational waves (Fig. 1b).
We predict the properties and detectability of the breakout emission in this
model, applied to possible electromagnetic counterparts of gravitational wave
events.

Figure 3.—  The SED for emission from typical merging black holes in AGN
disks. Same as Fig. 5, but the results assume common black hole mergers
($m=60~{}{\rm M}_{\odot}$ and $d_{\rm L}=500~{}{\rm Mpc}$) at
$R=1~{}{\rm{pc}}$ with $f_{\rm acc}=15$ (panel a) and at
$R=10^{-2}~{}{\rm{pc}}$ with $f_{\rm acc}=1$ (panel b). Dotted orange line
shows the sensitivity for Fermi GBM. The integration time is limited by
$t_{\rm duration}$ presented in the legend, while $t_{\rm
int}=10^{4}~{}{\rm{s}}$ (panel a) and $30~{}{\rm{s}}$ (panel b) are adopted
for the sensitivities of Swift BAT.  Figure 4.—  The timescales and luminosity
of shock breakout emission produced around merging black holes as a function
of the distance from the supermassive black hole ($R$). (a): The delay
timescale ($t_{\rm delay}$, black) and the observed duration of breakout
emission ($t_{\rm duration}$, orange). (b): The breakout luminosity by thermal
emission ($L_{\rm breakout}$, black) and that by the non-thermal emission at
its peak frequency ($L_{\rm non-thermal}$, orange) and at $h\nu=2~{}{\rm{eV}}$
($L_{\rm non-thermal}(2~{}{\rm{eV}})$, red), including Doppler beaming. The
non-thermal emission at $h\nu=2~{}{\rm{eV}}$ is also shown for the models with
lower efficiencies of magnetic field amplification ($\epsilon_{\rm{B}}=0.01$;
thin solid blue) and electron acceleration ($\epsilon_{\rm{e}}=0.03$: thin
dashed blue). The black hole locations adopted in the fiducial model are
indicated with filled circles superposed on the lines. The dashed holizontal
lines present the observed values for $t_{\rm delay}$, $t_{\rm duration}$, and
$L$ at $\sim 2~{}{\rm{eV}}$.  Figure 5.—  The SEDs for emission concerning the
putative associations with gravitational wave events. The SED for non-thermal
(solid black) and thermal (solid brown) emission from a merging black hole at
$R=2.5~{}{\rm{pc}}$ (panel a) and $R=10^{-4}~{}{\rm{pc}}$ (panel b). Solid and
dashed blue lines represent emission from the host AGN and its variability
amplitude, respectively. The dotted cyan, dark green, gray, purple, red, light
green, pink, and orange lines mark the sensitivity of ZTF, the Vera Rubin
Observatory and the Roman space telescope, Swift BAT (with the integration
time of $t_{\rm int}=10^{5}~{}{\rm{s}}$), HiZ-GUNDAM ($t_{\rm
int}=10^{4}~{}{\rm{s}}$) and Einstein Probe ($t_{\rm int}=10^{3}~{}{\rm{s}}$),
Chandra ($t_{\rm int}=2\times 10^{4}~{}{\rm{s}}$) and XMM-Newton ($t_{\rm
int}=10^{4}~{}{\rm{s}}$), NuSTAR ($t_{\rm int}=10^{6}~{}{\rm{s}}$) and FORCE
($t_{\rm int}=10^{6}~{}{\rm{s}}$), WISE ($t_{\rm int}=10^{4}~{}{\rm{s}}$), and
JWST ($t_{\rm int}=10^{4}~{}{\rm{s}}$), respectively. The red points in panel
(a) correspond to the observed luminosity of ZTF19abanrhr in the $r$ and $g$
bands assuming $g-r\sim 0.15$, and those in panel (b) represent the observed
luminosity of GW150914-GBM with $1$–$\sigma$ error bars. The dotted orange
curves mark the range of sensitivity of INTEGRAL SPI-ACS assuming $d_{\rm
L}=410~{}{\rm Mpc}$ (Savchenko et al., 2016). For thermal emission in panel
(b), we use the spectral shape given as the purple line in Fig. 4 in Ito et
al. (2020a) boosted to a higher energy by $\gamma_{\rm sf,f}=5$.

###  3.2. Properties of breakout emission

We first predict the characteristic properties of the breakout emission from a
merging black hole. To understand the detectability of the breakout emission
from common black hole mergers, Fig. 3 shows the parameter dependence of the
SED of the breakout emission adopting $m=60~{}{\rm M}_{\odot}$, $d_{\rm
L}=500~{}{\rm Mpc}$, and $R=1~{}{\rm{pc}}$ (panel a) and
$R=10^{-2}~{}{\rm{pc}}$ (panel b). We assume $f_{\rm acc}=15$ for
$R=1~{}{\rm{pc}}$ and $f_{\rm acc}=1$ for $R=10^{-2}~{}{\rm{pc}}$, since the
accretion rate can be enhanced by the recoil kick if the delay timescale is
long enough that shocked gas can accrete onto a black hole before the breakout
($\S$B). For mergers at $R=1~{}{\rm{pc}}$, non-thermal emission can be
discovered by ZTF, the Vera Rubin Observatory, the Roman space telescope, HiZ-
GUNDAM, and the Einstein Probe, and also be detected by follow-up observations
by Chandra, XMM-Newton, NuSTAR, FORCE, WISE, and JWST once the source
direction is localized by ZTF. Additionally, thermal emission from
$R=1~{}{\rm{pc}}$ and $R=10^{-2}~{}{\rm{pc}}$ can be detected by the Swift
Burst Alert Telescope (BAT). However, for $R=10^{-2}~{}{\rm{pc}}$, the
detection probability is reduced by the beaming effect by a factor of $\sim
1/\gamma_{\rm sf,f}^{2}\sim 1/16$ compared to cases with $\gamma_{\rm
sf,f}\sim 1$, since the opening angle of the jet at breakout of
$\theta_{j}\sim 0.08(<\theta_{0})$ (Bromberg et al., 2011) is smaller than
$1/\gamma_{\rm sf,f}$, where $\gamma_{\rm sf,f}$ is the final Lorentz factor
of the shocked gas. In addition, Swift BAT needs to direct to the source
beforehand since there is not enough time for the telescope to redirect due to
the short delay time ($\sim 10^{2}~{}{\rm{s}}$) and duration ($\sim
30~{}{\rm{s}}$), which further reduces the observational probability by a
factor of $\sim 10$. Thus, for mergers in $R\gtrsim{\rm{pc}}$, optical surveys
and/or future X-ray monitors can detect electromagnetic counterparts, as
suggested for ZTF19abanrhr, and then, infrared and X-ray pointing facilities
can detect it later. For $R\lesssim 10^{-2}~{}{\rm{pc}}$, electromagnetic
counterparts are difficult to discover unless they are bright enough to be
detected by the Fermi GBM and/or the INTEGRAL SPI-ACS as suggested for
GW150914-GBM and LVT151012-GBM.

In Fig. 3 a, the luminosity from the host AGN in the relevant energy range is
given by

$\displaystyle\nu L_{\rm AGN}(\nu)\sim 10^{44}~{}{\rm erg/s}(M/10^{8}~{}{\rm
M}_{\odot})({\dot{M}}c^{2}/L_{\rm Edd})(f_{\rm cor}/10)^{-1}\,.$ (38)

where $f_{\rm cor}$ is the correction fraction from the total luminosity to
the luminosity at some frequency. As depicted by the solid blue lines in Figs.
3 and 4, we assume that $f_{\rm cor}\sim 5$ at $c/\nu=4400~{}$Å and
extrapolate the luminosity for $10^{12}~{}{\rm Hz}\lesssim\nu\lesssim
10^{15}~{}{\rm Hz}$ using the cyan or blue points in Fig. 7 of Ho (2008)
depending on the assumed Eddington ratio, and $f_{\rm cor}\sim 10$ in
$0.1~{}{\rm keV}\leq h\nu$ (Ho, 2008; Trakhtenbrot et al., 2017; Duras et al.,
2020) with the upper exponential cut off at $300~{}{\rm keV}$ (Ricci et al.,
2018). We also assume that the fraction of the variable luminosity compared to
the average luminosity in optical bands with $t_{\rm duration}\lesssim
0.1~{}{\rm yr}$ is $f_{\rm var}\lesssim 0.1$ (Kozłowski, 2016) and that in
X-ray bands is $f_{\rm var}\sim 0.3$ (Soldi et al. 2014; Maughan & Reiprich
2019, dashed blue lines).

###  3.3. ZTF19abanrhr and GW190521

####  3.3.1 Observed properties

We overview observed properties for ZTF19abanrhr. Graham et al. (2020)
reported that the optical transient, ZTF19abanrhr, from the AGN ${\rm
J124942.3+344929}$ at $z=0.438$ is possibly associated with GW190521. The
optical flare in ZTF19abanrhr began to exceed the persistent flux from the AGN
activity after $\sim 18$ days (in the rest frame) from the occurrence of the
merger event, GW190521. The peak luminosity of ZTF19abanrhr is $\sim 8\times
10^{44}\,{\rm erg/s}$ in the $g$ and $r$ bands with the observed duration of
$\sim 28\,{\rm days}$ (in the rest frame). Before and after the peak, the
luminosity increases and decreases roughly exponentially. The change in the
slope of the luminosity is somewhat shallower in the latter. The color is
consistent with a constant with $g-r\sim 0.2$, while it appears to be reduced
to $g-r\sim 0.15$ during the flare. The mass of the supermassive black hole
hosting the AGN ${\rm J124942.3+344929}$ is $M\sim 10^{8}$–$10^{9}\,{\rm
M}_{\odot}$, and the luminosity of the AGN in units of the Eddington
luminosity is $\sim 0.02$–$0.2$ (Graham et al., 2020). The mass of the merged
remnant is $m\sim 150\,{\rm M}_{\odot}$ and its dimensionless spin is $a_{\rm
BH,rem}\sim 0.7$ (Abbott et al., 2020a). The inclination angle (angle between
the total angular momentum of the merged binary with respect to the line-of-
sight) is estimated to be $\sim 40$–$60^{\circ}$ (Abbott et al., 2020a;
Calderón Bustillo et al., 2021), disfavoring a transient with high luminosity
due to a significant beaming of emission from the Blandford-Znajek jet or its
shock, since the spin of a merged remnant, and accordingly the Blandford-
Znajek jet, are expected to be aligned with the total angular momentum of the
merged binary. Note that the inclination angle estimated by Gayathri et al.
(2022) is different from that above. Thus, the possible beaming of emission
associated with gravitational wave events may be useful to constrain the
parameters of black hole mergers observed by gravitational waves.

The association significance of ZTF19abanrhr to GW190521 is under debate
(Ashton et al., 2021; Palmese et al., 2021; Calderón Bustillo et al., 2021).
Ashton et al. (2021) analyzed that a chance occurrence probability of the
association of GW190521 and ZTF19abanrhr is between $\sim 8\%$–$50\%$. Palmese
et al. (2021) estimated that the chance occurrence probability is $\sim 4\%$
and $\sim 70\%$ by using a damped random walk model adopted in Graham et al.
(2020) and a more general structure function, respectively. Here, a damped
random walk model well describes structures of flares from AGN disks for the
duration of flares to be less than $\sim{\rm yr}$ (Kozłowski, 2016). Calderón
Bustillo et al. (2021) found that with different priors (uniform in inverse of
a mass ratio, $1$–$4$, instead of uniform in a mass ratio, $1/4$–$1$), the
chance occurrence probability is reduced to $\sim 1\%$, due to the
modification for the position of GW190521 in different preferences for the
mass ratio. Currently, the association is considered to be tentative. There
are also several independent suggestions that GW190521 originated in an AGN
disk (Yang et al., 2020; Tagawa et al., 2021a), especially if the suggested
high eccentricity is true (Tagawa et al., 2021b; Samsing et al., 2022; Romero-
Shaw et al., 2020; Gayathri et al., 2022).

####  3.3.2 Model for ZTF19abanrhr

We here discuss whether the possible association of ZTF19abanrhr to GW190521
can be explained by emission from black holes merging in an AGN disk. We focus
on this association here, but note that other optical associations were
reported (Graham et al., 2022) after the present manuscript was being prepared
for submission. These additional sources will be discussed in a separate
follow-up paper. We adopt parameter values to match the observed properties of
this source ($\S$2.4). To propose a consistent model for this association, it
is crucial to simultaneously explain (I) all of the several observed
properties of the flare and (II) the reason why there is no bright emission
before the merger.

(I) Several properties of the flare are obtained by the ZTF observations,
which include its luminosity in the optical bands, its color, and its time
evolution ($\S$3.3.1). Fig. 4 shows the properties of thermal and non-thermal
emission at the shock breakout of a jet produced from merging black holes. The
discontinuities at distances from the central supermassive black hole of
$R\sim 1~{}{\rm{pc}}$ and $\sim 0.01~{}{\rm{pc}}$ are, respectively, due to
the transition from sub-relativistic to relativistic regimes and to gap
formation around black holes by their gravitational torques. In this model,
the delay timescale ($t_{\rm delay}\sim 19~{}{\rm days}$, Fig. 4a) and the
observed duration of the breakout emission for non-thermal photons ($t_{\rm
duration}\sim 23~{}{\rm days}$, panel a) are, respectively, comparable to the
observed delay time and the duration (dotted gray lines), provided the merger
occurs at $R\sim 2.5$ pc. In the breakout emission, the delay timescale is
calculated by the time for the shock to reach the edge of the AGN disk (Eq.
2.2), while the duration is calculated by the time that photons diffuse out
from a breakout shell to the edge of the AGN disk (Eq. 2.2, see $\S$2.2).
Since the observable properties especially the timescales are mostly
influenced by $R$ (Fig. 4), flares similar to previous one (ZTF19abanrhr if it
was real) are expected to be observed in the future if there is a hot spot for
the merging location ($\S$4). Fig. 5 shows the spectral energy distribution
(SED) of the emission at $R=2.5~{}{\rm{pc}}$. The luminosity at $\sim 6\times
10^{14}~{}{\rm Hz}$, arising from synchrotron emission from the non-thermal
electrons accelerated at the forward shock, is roughly consistent with the
observed value of $\nu L_{\nu}\sim 8\times 10^{44}~{}{\rm erg/s}$. The power-
law spectral slope ($s$, $\nu L_{\nu}\propto\nu^{s}$) of the synchrotron
emission is given by $s=-0.25$, while the observed slope during the flare
corresponding to the color $g-r=0.15$–$0.2$ is $s=-(0.45$–$0.6)$ (Graham et
al., 2020). Since the contribution of the flare to the background AGN
luminosity is $\sim 30\%$ and the slope of the background emission is
$s\sim-0.6$ ($g-r=0.2$), the slope for the combination of the flare and the
background emission is estimated to be $\sim-0.5$ at the observed wavelengths
of $~{}470~{}{\rm nm}$ and $\sim 650~{}{\rm nm}$. Hence the synchrotron
emission is roughly consistent with the observed color during the flare. The
remarkable point of the estimates is that the delay time, the duration, and
the luminosity of ZTF19abanrhr are well reproduced for the same value of the
location within the disk, that is $R=2.5~{}{\rm{pc}}$ (Fig. 4), although the
luminosity is somewhat adjusted ($\S$2.4). For the non-thermal emission, the
peak of the lightcurve likely comes when the forward shock reaches the edge of
the AGN disk (Perna et al., 2021a). A more quantitative estimate of the
lightcurves for the non-thermal emission will require numerical simulations,
which is beyond the scope of this paper (see $\S$2.3).

(II) Another issue, which needs to be resolved to claim the association, is
that intense emission should be launched only after the merger of GW190521. In
the models for emission from the shock around the jet, this requires (a) the
existence of cold gas in the direction of the jet propagation from the merger
remnant and (b) gas accretion within the delay between GW190521 and
ZTF19abanrhr. The fact that the emission is produced only after the merger has
not been addressed or fully explained in previous models for the association.
Recoil kicks might offer a solution. Graham et al. (2020) proposed that the
merger remnant moves to unperturbed dense gas as a result of the recoil kick
and begins to accrete. However, it is unclear why circumbinary gas should not
be present and power bright emission already prior to the merger (Farris et
al., 2015; Tang et al., 2018). One possibility is that a cavity pre-exists
around the remnant, carved out by radiative and/or mechanical feedback from
the merging black holes; the merger remnant then moves to unshocked regions at
the boundary of this cavity by recoil kicks (Kimura et al., 2021). On the
other hand, with the parameters ($R\gtrsim{\rm{pc}}$) constrained above, it
takes too long ($\gtrsim{\rm yr}$) to reach replenished gas, compared to the
$\sim$ month delay between GW190521 and ZTF19abanrhr ($\S$4.1). As discussed
above, the model for breakout emission can explain why bright emission is
observed only after the mergers (Fig. 1b). We conclude that breakout emission
from shocks driven by a Blandford-Znajek jet can explain the properties of the
transient, ZTF19abanrhr, possibly associated with GW190521.

We further predict how the breakout emission similar to ZTF19abanrhr can be
observed by various telescopes in the future. We employ parameters to explain
ZTF19abanrhr as a fiducial model (see $\S$2.4). Fig. 5 shows the SED for the
breakout emission assuming a luminosity distance $d_{\rm L}\sim 3~{}{\rm Gpc}$
to the event together with the sensitivity curves of various equipments. Non-
thermal emission can be detected by ZTF, the Vera Rubin Observatory, and the
Roman space telescope (solid black, dotted cyan, and dotted green lines in
Fig. 5; see $\S$C for their properties). Thermal emission can be detected by
future wide-field X-ray surveys, HiZ-GUNDAM and the Einstein Probe if breakout
emission is produced closer to us (e.g. $d_{\rm L}\lesssim{\rm Gpc}$).
Chandra, XMM-Newton, Nuclear Spectroscopic Telescope Array (NuSTAR), Focusing
On Relativistic universe and Cosmic Evolution (FORCE, a future hard X-ray
telescope), Wide-field Infrared Survey Explorer (WISE), and the James Webb
Space Telescope (JWST) can detect emission from such merger remnants, if the
merger event is well localized by other observations, since the time spent for
directing to the source ($\sim{\rm day}$) is shorter than the duration of the
emission for mergers at $R\gtrsim{\rm{pc}}$ ($\sim 0.1~{}{\rm yr}$, Fig. 4b).
As the duration of the emission is the same for all the wavelengths, infrared
and X-ray bands, simultaneously with bright optical emission. We propose that
such follow-up observations in infrared and X-ray bands can be a smoking gun
signature of this scenario, leading to a first confirmation of the origin of
the black hole mergers. Then, we can further derive properties of the mergers’
environments as described for possible several associations above and in
$\S$2.

###  3.4. GW150914-GBM and LVT151012-GBM

####  3.4.1 Observed properties

We next describe the observed properties of GW150914-GBM. Associated with
GW150914, the Fermi GBM might have detected a transient of luminosity $\sim
2\times 10^{49}\,{\rm erg/s}$ at energies $\sim 10~{}{\rm keV}$–several ${\rm
MeV}$ with spectral shape $\nu L_{\nu}\propto\nu^{p_{\rm{s}}}$ with
$p_{\rm{s}}\sim 0.6$. The possible transient was observed $\sim 0.4\,{\rm{s}}$
after the gravitational wave event with duration of $\sim 1~{}$s. The signal-
to-noise ratio was estimated to be 5.1 with a false-alarm probability for the
association with GW150914 of 0.0022 (2.9$\sigma$) (Connaughton et al., 2016).
Several criticisms and/or issues were raised by a number of studies (Greiner
et al., 2016; Savchenko et al., 2016), and they are answered or discussed in
Connaughton et al. (2018). A possible problematic point is that the anti-
coincidence shield (ACS) of the Spectrometer on board INTEGRAL (SPI) put
constraints on the gamma-ray intensities in $75~{}{\rm keV}$–$2~{}{\rm MeV}$
in the direction of GW150914 at the merger (Savchenko et al., 2016). On the
other hand, since INTEGRAL is most sensitive at $\sim 100~{}{\rm keV}$ (Fig. 3
in Savchenko et al. (2016)), while the peak energy inferred for GW150914-GBM
is $\sim 3.5^{+2.1}_{-1.1}~{}{\rm MeV}$ and the estimated spectral index is
hard, it can be barely consistent with the no detection by INTEGRAL SPI-ACS
(Connaughton et al., 2018) (see Fig. 3b). Connaughton et al. (2018) discussed
that a weaker signal estimated by the analyses in Greiner et al. (2016) may be
preferred for the consistency with the no detection by INTEGRAL SPI-ACS.

Additionally, Bagoly et al. (2016) reported another possible electromagnetic
association with the gravitational wave event, LVT151012, with false alarm
probability of $0.037$. LVT151012 is a tentative gravitational wave event for
the merger of $23^{+18}_{-5}~{}{\rm M}_{\odot}$ and $13^{+4}_{-5}~{}{\rm
M}_{\odot}$ at $d_{L}=1100^{+500}_{-500}~{}{\rm Mpc}$ with a false alarm rate
of $0.44~{}{\rm yr}^{-1}$ and a false alarm probability of $0.02$ (Abbott et
al., 2016b). The burst looks similar to GW150914-GBM. It has a similar flux,
it occurred within $\sim~{}{\rm{s}}$ from the gravitational wave event with
duration of $\sim{\rm{s}}$, and the peak energy is found to be between
$130~{}{\rm keV}$ and $3.5~{}{\rm MeV}$.

####  3.4.2 Model for GBM events

We also find that the properties, including the luminosity, delay time,
duration and color of the transients suggested to be associated with GW150914
and LVT151012 can be explained by thermal emission from black hole mergers at
$R=10^{-4}~{}{\rm{pc}}$ and $R\lesssim 10^{-4}~{}{\rm{pc}}$, respectively.

A model for this candidate source needs to explain the observed hard spectrum
(red circles in Fig. 5b). A possible process reproducing the spectrum is
thermal emission. The peak energy of the transient is about $\sim 3~{}{\rm
MeV}$, while annihilation of gamma rays prohibits strong thermal emission
above $\sim{\rm MeV}$ in the rest frame (Budnik et al., 2010; Ito et al.,
2020a). To explain the strong emission at $\sim 3~{}{\rm MeV}$, the observer
needs to be within a beaming direction of the shock at breakout with a Doppler
shift by a factor of $\sim 4$–$10$ (Fig. 5b). The spectral shape of the
thermal emission in Fig. 5b is adopted from that calculated by Ito et al.
(2020a) considering $\gamma\gamma$ annihilation. For these events, we assume
the Bondi-Hoyle-Lyttleton rate for the accretion rate onto the merged remnant
($\S$2.1).

We constrain model parameters as follows. First, to produce the spectral peak
energy, we adopt $\gamma_{\rm sf,f}\sim 5$ reflecting the discussion above.
The observed delay timescale of $t_{\rm delay}\sim 0.4~{}{\rm{s}}$ and
$\gamma_{\rm sf,f}\sim 5$ together can be used to constrain the scale height
of the AGN disk to $H_{\rm AGN}\sim 3\times 10^{11}~{}{\rm cm}(\gamma_{\rm
sf,f}/5)^{2}(t_{\rm delay}/0.4~{}{\rm{s}})$ (see Eq. 7). The observed breakout
luminosity ($L_{\rm breakout}$) is proportional to the jet power ($L_{\rm j}$)
and the accretion rate onto the black hole ${\dot{m}}$ (Eq. 1), which is
adjusted by changing the AGN density at the location of the black hole
($\rho_{\rm AGN}$) and the supermassive black hole mass ($M$) given the black
hole mass ($m$), $R$, and $\gamma_{\rm sf,f}$ (Tagawa et al., 2022b).
Additionally, the initial opening angle of the jet $\theta_{0}$ at launch can
be determined to reproduce $\gamma_{\rm sf,f}$ with given $H_{\rm AGN}$,
$\rho_{\rm AGN}$, and $L_{\rm j}$ (Bromberg et al., 2011). By adopting
$\rho_{\rm AGN}=4\times 10^{-6}~{}{\rm g/cm^{3}}$, $M=10^{6}\,{\rm
M}_{\odot}$, $R=10^{-4}~{}{\rm{pc}}$, and $\gamma_{\rm sf,f}=5$, we can derive
$\theta_{0}=0.06$, $H_{\rm AGN}\sim 3\times 10^{11}$, $L_{\rm
j}=10^{48}~{}{\rm erg/s}$, $L_{\rm breakout}=1.5\times 10^{49}~{}{\rm erg/s}$
(including Doppler beaming), the delay time $t_{\rm delay}=0.4~{}{\rm{s}}$
(using Eq. 7), and the duration $t_{\rm duration}=0.2~{}{\rm{s}}$ (using Eq.
11). Note that $t_{\rm duration}=0.2~{}{\rm{s}}$ is shorter than the inferred
total duration of $\sim 1.0~{}{\rm{s}}$. On the other hand, the duration of
the phase with the above bright luminosity is $\lesssim 1/4~{}{\rm{s}}$ (Fig.
7 in Connaughton et al. (2016)), consistent with the value derived here (see
also Bagoly et al. (2016)). The light-curve is expected to grow super-
exponentially at breakout (Sapir et al., 2011), and after breakout emission,
the luminosity decays as a power-law in time (Nakar & Sari, 2012). The
inferred density ($\rho_{\rm AGN}=4\times 10^{-6}~{}{\rm g/cm^{3}}$) together
with $R=10^{-4}~{}{\rm{pc}}$, $H_{\rm AGN}\sim 3\times 10^{11}$, and
$M=10^{6}\,{\rm M}_{\odot}$ is realized for an AGN disk with a common
accretion rate of ${\dot{M}}\sim 0.2~{}{L}_{\rm Edd}/c^{2}$ and a viscous
parameter of $\alpha_{\rm AGN}=0.1$. Thus, even without identifying a host
AGN, if our model correctly predicts the electromagnetic emission, we can
derive properties of the AGN, such as $H_{\rm AGN}$, $\rho_{\rm AGN}$, and
$R$, where black hole mergers often occur. The spectral energy distribution
(SED) adjusted to explain GW150914-GBM is presented in Fig. 5b.

Although we consider the collision of the jet with the inner regions of the
AGN disk, it may also be possible to produce similar emission by considering
the collision of the jet with a circum-black hole disk, whose size is $\sim
3\times 10^{11}~{}{\rm cm}$ for the black hole mass of $m\sim 60~{}{\rm
M}_{\odot}$ and the accretion rate onto the black hole of ${\dot{m}}\sim
10^{4}~{}l_{\rm Edd}/c^{2}$ assuming a slim disk model (Abramowicz et al.,
1988).

The properties of the transient possibly associated with LVT151012 ($\S$3.4.1)
can also be explained by a similar parameter set, while the higher luminosity
(due to the larger distance by a factor of $\sim 2$) and the lower peak energy
(roughly around MeV) require a higher $\rho_{\rm AGN}$ and a lower
$\gamma_{\rm sf,f}$. These can be explained by a merger at smaller $R$, where
$\rho_{\rm AGN}$ is higher and a lower $\gamma_{\rm sf,f}$ is preferred to
reproduce the delay time ($\sim{\rm{s}}$) due to a smaller $H_{\rm AGN}$.

In summary, in the model of emission from merging black holes in an AGN disk,
there exists a set of parameters for which the properties of ZTF19abanrhr,
GW150914-GBM, and LVT151012-GBM can be explained.

##  4\. Discussions

In this section, we discuss the probability of observing breakout emission
from merging black holes, the possible obscuration by a dust torus, the
parameter dependence of the results, and the merging location within the AGN
disc. Specifications of several telescopes which might be able to detect the
breakout emission at different wavelengths are summarized in Appendix C.

###  4.1. Probability of observing breakout emission from merging black holes

We suggest that breakout emission from merging black holes may be easier to
discover compared to that from solitary black holes due to the following
reasons: (A) By using gravitational wave observations, constraints on the
spatial localization and the timing of flares are modestly and significantly
improved, respectively, thus removing the need to continuously monitor large
sky areas to search for rare flares. (B) The accretion rates onto merger
remnants may be enhanced due to recoil kicks (see below). (C) The spin
magnitudes of black holes are enhanced to $a_{\rm BH,rem}\sim 0.7$ by the
merger (Buonanno et al., 2008), increasing the jet luminosity. Although the
spin magnitudes of solitary black holes ($a_{\rm BH,iso}$) are highly
uncertain, spin-up by accretion and spin-down by the Blandford-Znajek jet may
be roughly equal at around $a_{\rm BH,iso}\lesssim 0.3$ (e.g., Fig. 10 of
Narayan et al. 2022). Assuming $a_{\rm BH,iso}\sim 0.3$, the jet luminosity is
enhanced by a factor of $\sim(a_{\rm BH,rem}/a_{\rm BH,iso})^{2}\sim 5(a_{\rm
BH,rem}/0.7)^{2}(a_{\rm BH,iso}/0.3)^{-2}$. (D) The merging black holes may be
more massive than isolated black holes as inferred from gravitational wave
observations (The LIGO Scientific Collaboration et al., 2021), producing
brighter emission due to higher accretion rates.

In order for emission to be produced shortly after mergers (roughly within a
timescale comparable to the duration of the flare) and be discovered using
gravitational wave observations, there are several key elements which are
needed. First, (a) gas needs to be present along the direction of propagation
of the jet (black hole spin) of the merged remnant, or otherwise thermal and
non-thermal emission from dissipation of the jet power would not be expected.
The shape of the ejected region is cylindrical with an aspect ratio of
$\sim\theta_{\rm 0}f_{\rm ext}\sim 1(\theta_{0}/0.3)(f_{\rm ext}/3)$ (Tagawa
et al., 2022b), where $\theta_{0}$ is the opening angle of the injected jet,
$f_{\rm ext}$ is the fraction of the radial extent of the shocked regions
advanced after the breakout from the AGN disk over that at the breakout, and
$f_{\rm ext}\sim 3$ is motivated for a spherical explosion in a background
with an exponential gas density profile (Olano, 2009). Hence, gas presumably
exists at the merger except in the cylindrical regions that have been ejected
before the merger in previous episodes of jets launched from the progenitor
black holes. Since the angular momentum directions of binaries are predicted
to be random at mergers in AGN disks due to frequent binary-single
interactions (Tagawa et al., 2020a; Samsing et al., 2022) and inhomogeneities
in the outer regions of the AGN disk, the directions of the spins of the
merged remnants and that of the jet propagation are plausibly randomized at
mergers. Thus, gas often exists in the direction of the spin of the merged
remnants (Fig. 1b), with a probability of $P_{\rm col}\sim{\rm cos}({\rm
atan}(\theta_{\rm 0}f_{\rm ext}))\sim 0.7$ for $\theta_{\rm 0}\sim 0.3$ and
$f_{\rm ext}\sim 3$. If we assume $\theta_{0}\sim 0.03$–$0.5$ (Berger, 2014;
Hada et al., 2018) and $f_{\rm ext}\sim 1$–$5$ (Olano, 2009), $P_{\rm col}$
ranges in the interval $\sim 0.4$–$1$.

The second issue is that (b) gas needs to be accreted onto the merged remnant
for a flare to occur in association with a black hole binary merger. Here, we
assume that the flare can be considered to be associated with a gravitational
wave event when the delay timescale is roughly comparable to or less than the
duration of the flare, since it is probably difficult to find an association
with a short duration and a long delay time. In Kimura et al. (2021), it is
proposed that electromagnetic counterparts may be produced by moving to
unshocked regions by the gravitational wave recoil kick after merger. Usually,
gas replenishment by the recoil kick takes a longer time compared to $t_{\rm
delay}$ and $t_{\rm duration}$. The kick velocity due to gravitational wave
radiation ($v_{\rm rk}$) is on the order of $\sim
10^{2}$–$10^{3}~{}{\rm{km/s}}$ and gas is depleted within $r_{\rm
dep}\sim~{}f_{\rm ext}H_{\rm AGN}\theta_{\rm 0}\sim 10^{16}\,{\rm
cm}(R/{\rm{pc}})$ (Tagawa et al., 2022b), where $H_{\rm AGN}$ is the scale
height of the AGN disk at $R$. Then, the kicked black hole crosses the cavity
on the timescale of

$\displaystyle t_{\rm arr}\sim r_{\rm dep}/v_{\rm rk}\sim 30\,{\rm
yr}~{}(v_{\rm rk}/3\times 10^{2}\,{\rm{km/s}})^{-1}(R/3~{}{\rm{pc}}).$ (39)

Such a long delay time ($t_{\rm arr}$) compared to the short duration (orange
line in Fig. 2 a) would not be regarded as association in observations. In
ZTF19abanrhr, the observed delay time of $\sim 0.05~{}{\rm yr}$ is much
shorter than the arrival timescale in Eq. (39). To expect a flare with a short
delay whose timescale is comparable to the observed delay time, the merging
black holes need to have circum-black hole disks at the merger. Adopting the
model in Tagawa et al. (2022b), we derive that $P_{\rm active}\sim 0.03$ at
$R\sim 3~{}{\rm{pc}}$, where $P_{\rm active}$ is the probability that a merged
remnant accompanies a circum-black hole disk and a jet, although the active
probability strongly depends on model parameters. On the other hand, if
circum-black hole disks play a significant role for driving mergers, merged
remnants may tend to accompany circum-black hole disks (Farris et al., 2015;
Tang et al., 2018). Bartos et al. (2017) assessed mergers to be significantly
facilitated by circum-black hole disks, while Tagawa et al. (2020b) found the
opposite by applying the updated theory of gap formation (Duffell et al.,
2014). Note that both studies assumed that the accretion onto black holes is
limited at around the Eddington rate, which may underestimate the effect on
binary evolution by circum-black hole disks because the disk mass can be much
higher for highly super-Eddington accretion disks; hence the influence of
torques from circum-black hole disks on mergers has not been quantitatively
understood. Since $P_{\rm active}$ is estimated to be larger than $\sim 0.01$
(Tagawa et al., 2022b), and significant facilitation of mergers by circum-
binary disks may enhance $P_{\rm active}$ up to 1, we assume as a rough
estimate that $P_{\rm active}$ ranges within $\sim 0.01$–$1$.

Third, (c) the breakout emission needs to be bright enough to overshine the
host AGN emission and (d) bright enough to be detected by current (and future)
facilities. However, the probabilities that the breakout emission is brighter
than the host AGN emission ($P_{\rm bright}$) and that can be detected by
facilities ($P_{\rm det}$) are uncertain and difficult to constrain
theoretically as the luminosity is influenced by several uncertain parameters
(such as $\epsilon_{\rm{B}}$, $\epsilon_{\rm{e}}$, $f_{\rm acc}$, and
$\rho_{\rm AGN}$) and also transients with a high shock Lorentz factor
($\gamma_{\rm sf,f}$), as suggested for GW150914-GBM, are rarer by a factor of
($\sim\gamma_{\rm sf,f}^{2}$) to be in the beaming direction of the shock.

Finally, (e) black hole mergers need to occur in AGN disks. Although this
probability ($P_{\rm AGN}$) has been highly debated (Gayathri et al., 2021b),
it has been less constrained and ranges from 0 to 1.

Considering these factors and adopting the fiducial values for $P_{\rm col}$
and $P_{\rm active}$ and $P_{\rm bright}P_{\rm det}P_{\rm AGN}=1$, the
probability for observing breakout emission associated with mergers in AGN
disks is $P_{\rm association}\sim~{}P_{\rm col}P_{\rm active}P_{\rm
bright}P_{\rm det}P_{\rm AGN}\sim 0.02(P_{\rm col}/0.7)(P_{\rm
active}/0.03)(P_{\rm bright}P_{\rm det}P_{\rm AGN}/1)$. According to this
estimate, if $1$–$9$ transients implied in Results are actually associated
with the mergers among O(90) mergers discovered by LIGO/Virgo/KAGRA, as
modeled above, $P_{\rm bright}P_{\rm det}P_{\rm AGN}$ can be constrained to be
high ($\gtrsim 0.5$). Note that there are significant uncertainties as $P_{\rm
col}P_{\rm active}$ ranges from $4\times 10^{-3}$ to $1$, and $P_{\rm
bright}$, $P_{\rm det}$, and $P_{\rm AGN}$ range from 0 to 1 as discussed
above. Nevertheless, we expect that future observations, such as by the Vera
Rubin Observatory and the Roman space telescope will significantly increase
$P_{\rm det}$ especially for dimmer breakout emission from less bright AGNs as
well as constrain $P_{\rm col}P_{\rm active}$, $P_{\rm bright}$, and $P_{\rm
AGN}$ due to their high sensitivity and wide field-of-view.

###  4.2. Possible obscuration by a dust torus

If abundant dust exists along the line of sight to the observer, emission in
the optical to soft X-ray bands can be absorbed. Due to the existence of a
dust torus around AGNs, the breakout emission in these bands may hence be
largely absorbed. In particular, our model predicts a merger at $R\sim
2.5~{}{\rm{pc}}$ for ZTF19abanrhr, which is beyond the dust sublimation radius
($\sim 0.1~{}{\rm{pc}}$ for the AGN luminosity of $\sim 10^{44}~{}{\rm erg/s}$
in the ultraviolet bands, Barvainis 1987), suggesting the existence of dust in
this region. On the other hand, for $R\sim$a few pc, dust in regions above and
below the AGN disk plane is possibly cleared out by AGN disk winds (Wada et
al., 2016) and/or supernovae feedback (Wada et al., 2009) (see also
observations, Stalevski et al. 2017). If this is the case, dust in the disk is
usually confined to the disk plane and can be evaporated locally by the cocoon
feedback, as also confirmed in gamma-ray bursts (Waxman & Draine, 2000). The
breakout emission can therefore emerge without obscuration by dust in this
case. Conversely, measured properties (e.g. normalization and reddening) of
the breakout emission from $R\sim$a few pc may be useful to understand the
geometry of dust tori. Also, if the breakout emission is produced from black
holes in type II AGNs unlike the host of ZTF19abanrhr, the emission in optical
and soft X-ray bands is likely obscured. Since the emission is unobscured in
infrared and hard X-ray bands even if the emission is produced in
geometrically thick dust tori, simultaneous observations with infrared/hard
X-ray and optical bands will be able to constrain the configuration of dust
tori. Thus, wide-field infrared/hard X-ray facilities, such as Swift BAT, will
play important roles to find the electromagnetic counterparts produced in dust
tori. Additionally, due to the obscuration, $P_{\rm obs}$ is reduced by the
covering fraction of the dust torus. Although the influence of dust is
uncertain, our model can avoid the issue of dust obscuration to explain
ZTF19abanrhr for the reasons discussed above.

Figure 6.—  Similar to Fig. 2 but results for several choices of the
parameters are shown. Black, orange, blue, and red lines, respectively,
present results for the fiducial model, and the models with ${\dot{M}}_{\rm
in}=2~{}{\rm M}_{\odot}/{\rm yr}$, $M=10^{7}~{}{\rm M}_{\odot}$ and
${\dot{M}}_{\rm in}=1~{}{\rm M}_{\odot}/{\rm yr}$, and $f_{\rm acc}=5$. (a)
The delay time. (b) The duration (solid) and the breakout timescale (dashed).
(c) The shock (solid) and breakout (dashed) luminosities. (d) The absorption
(solid), minimum (dashed), and maximum (dotted) Lorentz factors for non-
thermal electrons. (e) The breakout temperature. (f) The dimensionless forward
shock velocity ($\gamma_{\rm FS}\beta_{\rm FS}$).

###  4.3. Parameter dependence

Fig. 6 shows the parameter dependence of the properties of breakout emission,
as shown in Fig. 2 for the fiducial model. For low ${\dot{M}}_{\rm in}$
(orange and blue lines) or low $f_{\rm acc}$ (red lines), the accretion rate
onto the BH is rather low, which lowers $L_{\rm j}$, $L_{\rm sh}$, and $L_{\rm
breakout}$ (panel c). At large $R$, for low ${\dot{M}}_{\rm in}$ (orange
lines), the luminosity is higher than that of the fiducial model. This is
because for low ${\dot{M}}_{\rm in}$ the AGN disk becomes geometrically thin
in outer regions, in which the jet head is less decelerated before the
breakout because the jet sweeps up a lower amount of gas. This enables a
higher speed for the shock. Note that the AGN disk density in outer regimes,
where the Toomre parameter satisfies $\sim 1$, only depends on the orbital
frequency, and is unaffected by ${\dot{M}}_{\rm in}$ (see Eq. 3 of Thompson et
al. 2005). Although there is the exception as explained above, to produce
bright emission (high $L_{\rm sh}$), high ${\dot{M}}_{\rm in}$ and/or high
$f_{\rm acc}$ are indeed required.

The fiducial model parameters and their influence on properties of the
breakout emission are listed in Table 1. There are several sensitive
parameters, $R$, $f_{\rm acc}$, $\epsilon_{\rm{B}}$, and $\epsilon_{\rm{e}}$.
$L_{\rm break}$ is significantly influenced by all of the four parameters,
while $t_{\rm delay}$ and $t_{\rm duration}$ are mostly influenced by $R$.
According to our model, $R$ can be constrained from the timescales, while
$f_{\rm acc}$, $\epsilon_{\rm{B}}$, and $\epsilon_{\rm{e}}$ can be constrained
if non-thermal emission is observed by multiple wavelengths.

Table 1 Fiducial values of our model parameters and their influence on properties of the breakout emission. For ”sensitivity”, we made a rough classification of the parameters influencing the results ($L_{\rm breakout}$, $t_{\rm delay}$, $t_{\rm duration}$) by more than two orders of magnitude (”very sensitive”), one order-of-magnitude (”sensitive”), a factor of $\gtrsim 2$ (”moderate”), and a factor of $\lesssim 2$ (”insensitive”). We also add the comment ”constrained” for parameters which are well constrained by the observations of ZTF19abanrhr.  Parameter |  Fiducial value |  Sensitivity to results
---|---|---
Radial distance of the black hole from the supermassive black hole |  $R=2.5\,{\rm{pc}}$ |  very sensitive on the three properties
Ratio of the black hole accretion rate to the gas capture rate |  $f_{\rm acc}=15$ |  sensitive on $L_{\rm breakout}$ and moderate on $t_{\rm delay}$ and $t_{\rm duration}$
Conversion efficiency of shock kinetic energy to magnetic/electron energy |  $\epsilon_{\rm{B}}=0.1$, $\epsilon_{\rm{e}}=0.3$ |  very sensitive on reducing $L_{\rm non-thermal}$
Jet energy conversion efficiency |  $\eta_{\rm j}=0.5$ |  sensitive on $L_{\rm breakout}$ and moderate on $t_{\rm delay}$ and $t_{\rm duration}$ but constrained
Mass of a stellar-mass black hole |  $m=150\,{{\rm M}_{\odot}}$ |  sensitive on $L_{\rm breakout}$ and moderate on $t_{\rm delay}$ and $t_{\rm duration}$ but constrained
Gas inflow rate from the outer boundary of the AGN disk |  ${\dot{M}}_{\rm in}=50~{}{L}_{\rm Edd}/c^{2}$ |  sensitive on the three properties but constrained
Mass of the supermassive black hole |  $M=10^{8}\,{{\rm M}_{\odot}}$ |  sensitive on the three properties but constrained
Opening angle at the base of the jet |  $\theta_{\rm 0}=0.3$ |  moderate on the three parameters
Angular momentum transfer parameter in outer regions of the AGN disk |  $m_{\rm AM}=0.5$ |  moderate on the three parameters
Viscous parameter of the AGN disk |  $\alpha_{\rm AGN}=0.1$ |  insensitive on the three parameters

###  4.4. Merging location

We have proposed that black hole mergers at $R\gtrsim{\rm{pc}}$, $R\sim
10^{-4}\,{\rm{pc}}$, and $R\lesssim 10^{-4}~{}{\rm{pc}}$ can explain the
properties of the transients possibly associated with GW190521, GW150914, and
LVT151012, respectively. Here we discuss whether such diversity in the
positions of mergers in AGN disks is expected. Mergers have been predicted to
be concentrated in gap-forming regions (Tagawa et al., 2020b, 2021a) and
migration traps (McKernan et al., 2020; Secunda et al., 2019; Yang et al.,
2020) in the inner AGN disk, or in the outer regions (Perna et al., 2021b) of
the disk (but may also occur throughout the disk, McKernan et al. 2020), which
have similar influence on the evolution and mergers of black hole binaries.
Gaps are predicted to form for massive black holes or in a thin disk (Tagawa
et al., 2020b), while migration traps may form at $\lesssim 10^{3}~{}R_{\rm
g}\sim 10^{-4}~{}{\rm{pc}}~{}(M/2\times 10^{6}~{}{\rm M}_{\odot})$ (Bellovary
et al., 2016) (but see Pan & Yang 2021), where $R_{\rm g}=GM/c^{2}$, $G$ is
the gravitational constant, and $c$ is the speed of light. If heating torques
are taken into account for the dynamics of gas around black holes (Hankla et
al., 2020), they can reverse the direction of migration at various distances
from the supermassive black hole, while the evolution may be further
complicated due to shocked bubbles (Tagawa et al., 2022b). Also, the migration
timescale without gaps is longer in outer regions especially around massive
supermassive black holes (Pan & Yang, 2021), in which mergers are also
expected at $R\gtrsim{\rm{pc}}$ (Perna et al., 2021b). As ZTF19abanrhr is
believed to be associated with an AGN powered by a supermassive black hole of
$M\sim 10^{8}\,{\rm M}_{\odot}$, inefficient migration (possibly due to long
migration timescale, gaps, or heating torques) may be consistent with mergers
at $R\gtrsim{\rm{pc}}$. Also, the mergers at $R\lesssim 10^{-4}\,{\rm{pc}}$,
suggested for the transients possibly associated with GW150914 and LVT151012,
may be consistent with mergers in migration traps of AGN disks around less
massive supermassive black holes with $M\sim 10^{5}$–$10^{6}\,{\rm
M}_{\odot}$. Conversely, if these transients are actually associated with
merging black holes in AGN disks, this places constraints on the positions of
gap-forming regions and/or migration traps, and yields information on the AGN
disks from the properties of the transients.

Since transients are reported for the massive mergers (including GW190521) or
the early events (GW150914 and LVT151012), the detection of transients may
require mergers of massive black holes and/or mergers at close distances from
us. Indeed, the breakout luminosity ($L_{\rm breakout}$), the jet power
($L_{\rm j}$), and the accretion rate onto a black hole (${\dot{m}}$) depend
on the mass of a black hole ($m$) as $L_{\rm breakout}\propto~{}L_{\rm
j}\propto{\dot{m}}\propto m^{2/3}$ (Eq. 1 of Tagawa et al. (2022b)), and the
flux is higher for closer events. Also, for transients with a longer delay
($t_{\rm delay}$), such as ZTF19abanrhr, enhancement of the accretion rate due
to recoil kicks ($\S$B) can be expected for massive remnants since the circum-
black hole disk mass within some radius $r$ from the black hole is $m_{\rm
CBD}(<r)\propto~{}m^{11/5}$ (Haiman et al., 2009).

##  5\. Summary

In this paper, we investigate properties of the EM emission associated with
black hole mergers in AGN disks. Both thermal and non-thermal radiation is
produced by shocks, through collisions between the AGN disc gas and a
Blandford-Znajek jet powered by the merger-remnant BH, in a direction
reoriented at the merger. We suggest that the emission emerging as the jets
break out from the optically thick AGN disc can self-consistently account for
several features of the tentative EM counterparts for the three GW events
ZTF19abanrhr, GW150914-GBM, and LVT151012-GBM. In particular, the model
accounts for the luminosity in the optical bands, as well as the timing and
duration of the emission. Our main results are summarized as follows.

  1. 1.

For mergers at a few pc from the central SMBH, the properties of ZTF19abanrhr
are reproduced by non-thermal emission. In this event, the accretion rate onto
the black hole is predicted to be enhanced by shocks due to recoil kicks.

  2. 2.

For mergers in closer vicinity of the central SMBH, the properties of
GW150914-GBM and LVT151012 are reproduced by thermal emission.

  3. 3.

The variety of distances from the central SMBH at which mergers are predicted
to take place is consistent with predictions of models for the AGN disk-
embedded BH population, in which the BHs migrate slowly and produce frequent
mergers.

To confirm this scenario, discoveries of electromagnetic counterparts for
future GW events are highly desired. Then, a smoking-gun signature, i.e.
bright infrared and X-ray emission concurrent with the optical transient, can
be detected by WISE, JWST, Chandra, XMM-Newton, NuSTAR, and FORCE, once host
galaxies are localized by optical observations, and/or Swift BAT could be able
to detect the breakout emission. If the scenario is confirmed, these systems
will be highly useful to improve our understanding of the evolution of compact
objects in AGN disks, the structure of AGN disks, plasma physics, and the
expansion history of the universe by helping to constrain the Hubble constant,
in addition to unveiling the origin of black hole mergers.

This work was financially supported by Japan Society for the Promotion of
Science (JSPS) KAKENHI grant Number JP21J00794 (HT) and 22K14028 (S.S.K.).
S.S.K. acknowledges the support by 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 NNX15AB19G and NSF grants AST-2006176 and AST-1715661. 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.
In Appendix, we describe the notations, the enhancement of the accretion rate
due to recoil kicks, and the properties of telescopes.

##  Appendix A Notation

The notations of variables are listed in Table 2.

Table 2 Notation.  Symbol |  Description |  Symbol |  Description
---|---|---|---
$r$ |  The distance from the stellar-mass black hole |  $H_{\rm AGN}$ |  The scale height of the AGN disk
$r_{\rm{b}}$, $r_{\rm ub}$ |  The radius within which all gas in a circum-black hole disk is bound and that beyond which all gas is unbound after the recoil kick |  $r_{\rm acc}$ |  The radius at which gas accretes onto the stellar-mass black hole on the breakout timescale
$R_{\rm g}$ |  The gravitational radius of a supermassive black hole $GM/c^{2}$ |  $r_{\rm k}$ |  The gravitational radius after a recoil kick $Gm/v_{\rm rk}^{2}$
$d_{\rm L}$, $d_{\rm L,3Gpc}$ |  The luminosity distance to the source and that in unit of $3~{}{\rm Gpc}$ |  $\Delta_{\rm shell}(\gamma)$ |  The shell width of electrons with $\gamma$ emitting non-thermal photons
$d_{\rm edge}$, $d_{\rm edge,BO}$ |  The distance to the surface of the AGN disk and that for the breakout emission |  $R_{\rm out}$ |  The outer boundary of the AGN disk
$m$ |  The mass of a stellar-mass black hole |  $m_{\rm CBD}(r)$ |  The mass of the circum-black hole disk within $r$ from a stellar-mass black hole
${\dot{m}}$, ${\dot{M}}$ |  The accretion rate onto a stellar-mass and a supermassive black hole, respectively |  ${\rm{m}}_{\rm BHL}$ |  The Bondi-Hoyle-Lyttleton rate
${\dot{m}}_{\rm Edd}$, ${\dot{M}}_{\rm Edd}$ |  The Eddington accretion rate onto a stellar-mass black hole and a supermassive black hole |  $l_{\rm Edd}$, $L_{\rm Edd}$ |  The Eddington luminosity of a stellar-mass and a supermassive black hole
$L$, $\nu L_{\nu}$ |  The luminosity and the luminosity at the frequency $\nu$ |  $L_{\rm j}$ |  The kinetic luminosity of the jet
$L_{\rm sh}$ |  The shock luminosity of the jet |  $L_{\rm breakout}$, $L_{\rm non-thermal}$ |  The breakout luminosity by thermal emission and that by the non-thermal emission at its peak frequency, respectively
${\tilde{L}}$ |  The ratio between the energy density of the jet and the rest-mass energy density of the surrounding medium |  $L_{\rm syn}$ |  The luminosity of synchrotron emission
$\eta_{\rm j}$ |  The conversion efficiency of mass to the jet |  $\eta_{\rm rad}$ |  The conversion efficiency of mass to radiation in an accretion disk
$Y_{\rm SSC}$, $Y_{\rm 2ndIC}$ |  The powers of synchrotron-self Compton and second order inverse Compton scattering compared to that of synchrotron emission |  $E_{\gamma,\rm BO}$ |  The energy of a photon
$t_{\rm duration}$, $t_{\rm delay}$ |  The duration and delay time of a transient |  $t_{\rm break}$,$t_{\rm break,rel}$, $t_{\rm GW,break}$ |  The breakout timescale of a non-relativistic and relativistic shock and gravitational wave, respectively
$t_{\rm diff}$ |  The diffusion timescale of photons out of the disk |  $t_{\rm ang}$ |  The angular timescale
$t_{\rm dyn,CBD}(r)$ |  The dynamical timescale of a circum-black hole disk at $r$ |  $t_{\rm vis}(r)$, $t_{\rm vis,sh}(r)$ |  The viscous timescale of a circum-black hole disk at $r$ before and after a shock arises due to a recoil kick
$t_{\rm int}$ |  The integration time of observations |  $t_{\rm arr}$ |  The timescale for a black hole arriving outside a cavity
$t_{\rm sph}$ |  The time of transition between the planar and spherical geometries of the breakout shell |  $t_{\rm dyn}$ |  The dynamical timescale of a shell
$t_{\rm trans}$ |  The transparent timescale |  $T_{\rm breakout}$ |  The breakout temperature
$\gamma_{\rm FS}$, $\gamma_{\rm sf}$, $\beta_{\rm FS}$, $\beta_{\rm sf}$ |  The Lorentz factor of the forward shock and shocked fluid, their dimensionless velocity |  $\gamma_{\rm FS,f}$, $\gamma_{\rm sf,f}$, $\beta_{\rm FS,f}$, $\beta_{\rm sf,f}$ |  The final Lorentz factor of the forward shock and shocked fluid, their dimensionless velocity
$\gamma_{\rm{e}}$ |  The Lorentz factor of a electron |  $\gamma_{\rm{m}}$, $\gamma_{\rm{max}}$, $\gamma_{\rm{a}}$, $\gamma_{\rm c}$ |  The minimum, maximum, absorption, and cooling Lorentz factors
Table 2continued  Symbol |  Description |  Symbol |  Description
---|---|---|---
$\nu$ |  The frequency of a photon |  $\nu_{\rm{syn}}$, $\nu_{\rm{m}}$, $\nu_{\rm{max}}$, $\nu_{\rm{a}}$,$\nu_{\rm{m,SSC}}$, $\nu_{\rm{a,SSC}}$ |  The synchrotron frequency for electrons with $\gamma_{\rm{e}}=1$, $\gamma_{\rm{m}}$, $\gamma_{\rm{max}}$, and $\gamma_{\rm{a}}$, and synchrotron self-Compton frequency for $\gamma_{\rm{m}}$ and $\gamma_{\rm{a}}$
$a_{\rm BH}$ |  The dimensionless spin parameter of a black hole |  $a_{\rm BH,rem}$, $a_{\rm BH,iso}$ |  The dimensionless spin parameter of a merged remnant and an isolated black hole
$f_{\rm{acc}}={\dot{m}}/{\dot{m}}_{\rm BHL}$ |  The enhancement factor of the accretion rate of the captured gas |  $f_{\rm ext}$ |  The fraction that the cocoon proceeds to the $r$-direction after the breakout
$P_{\rm col}$ |  The probability that a jet launched from a merged remnant collide with unshocked gas |  $P_{\rm active}$ |  The probability that a merged remnant is accreting immediately after merger
$P_{\rm obs}$ |  The probability that the emission from the breakout of the shock from the AGN disk is observable |  $v_{\rm rk}$ |  The recoil velocity due to anisotropic emission of gravitational waves at merger
$f_{\rm{b}}$, $f_{\rm ub}$ |  The factors for recoil kicks $f_{\rm{b}}=r_{\rm{b}}/r_{\rm k}$, $f_{\rm ub}=r_{\rm ub}/r_{\rm k}$ The fraction that the breakout emission is observable |  $f_{\rm accum}$ |  The fraction of gas within $r_{\rm ub}$ accumulated within $r_{\rm{b}}$ after the recoil kick
$f_{\rm inc}$ |  The factor that the gas mass within $r_{\rm{b}}$ is enhanced by the recoil shock |  $f_{\rm beaming}$ |  The beaming factor taking into account an angular effect
$f_{\rm corr}$ |  The correction factor from the total luminosity to the luminosity at some frequency |  $f_{\rm var}$ |  The ratio of the variable luminosity compared to the average luminosity at some frequency
$P_{\rm association}$ |  The probability that black hole mergers in AGN disks accompany electromagnetic observations |  $i$ |  The inclination angle between the jet and the orbital angular momentum of the AGN disk
$\theta_{\rm j}$, $\theta_{\rm 0}$ |  The jet opening angle and that at injection |  $\theta_{\rm rk}$ |  The angle of the recoil kick with respect to the plane of the circum-black hole disk
$\theta_{\rm c}$ |  The opening angle of the cocoon |  $\tau$ |  The optical depth of the AGN disk
$\kappa_{\rm CBD}$, $\kappa_{\rm AGN}$ |  The opacity of the circum-black hole disk and the AGN disk |  $\alpha_{\rm CBD}$ |  The viscous parameter for the circum-black hole disk
$\rho_{\rm AGN}$, $n_{\rm AGN}$ |  The density and the number density of the AGN disk |  $\rho_{\rm sf}$ |  The density of the shocked material
$e_{\rm sf}$ |  The energy density of the shocked material |  $B_{\rm sf}$ |  The magnetic field of the shocked material
$\xi$ |  The parameter representing the ratio of the mean free path to the Larmor radius of electrons |  $\tau_{q}$, $q$, $\gamma_{1}$, $\gamma_{2}$, $n_{\rm norm}$, $C_{q}$ |  The parameters for the optical depth by synchrotron self-absorption
$\theta_{\rm PA}$ |  The pitch angle between the magnetic field and the velocity of electrons |  $p$ |  The power-law slope for injected electrons accelerated by the first order Fermi process
$N(\gamma)d\gamma$ |  The number of electrons at $\gamma_{e}=\gamma\pm d\gamma/2$ |  $p_{\rm dm}$ |  The power-law slope for the mass profile of the circum-black hole disk
$s$ |  The power law slope for electromagnetic radiation at some frequency |  $e$ |  The elementary charge
$c$ |  The speed of light |  $G$ |  Gravitational constant
$k_{\rm{B}}$ |  The Boltzmann constant |  $h$ |  The Planck constant
$m_{p}$ |  The proton mass |  $a$ |  The radiation constant

##  Appendix B Enhancement of the accretion rate onto merger remnants

For accretion onto merged remnants, the accretion rate can be significantly
affected by the recoil kicks due to anisotropic radiation of gravitational
wave. We follow Rossi et al. (2010) for estimating the enhancement of the
accretion rate below. After the recoil kick, all gas in an accretion disk
keeps bound within the radius of

$\displaystyle r_{\rm{b}}$ $\displaystyle=f_{\rm{b}}r_{\rm
k}=\frac{f_{\rm{b}}Gm}{v_{\rm rk}^{2}}$ $\displaystyle=5\times 10^{12}~{}{\rm
cm}\left(\frac{f_{\rm{b}}}{0.2}\right)\left(\frac{m}{150~{}{\rm
M}_{\odot}}\right)\left(\frac{v_{\rm rk}}{300~{}{\rm{km/s}}}\right)^{-2},$
(B1)

while all gas is unbound beyond the radius of

$\displaystyle r_{\rm ub}=f_{\rm ub}r_{\rm k},$ (B2)

where $r_{\rm k}\equiv Gm/v_{\rm rk}^{2}$, $f_{\rm{b}}=[-{\rm cos}\theta_{\rm
rk}+({\rm cos}\theta_{\rm rk}^{2}+1)^{1/2}]^{2}$ and $f_{\rm ub}=[{\rm
cos}\theta_{\rm rk}+({\rm cos}\theta_{\rm rk}^{2}+1)^{1/2}]^{2}$ are the
factors related to the the kick direction (e.g. $f_{\rm{b}}\sim 0.2$ and
$f_{\rm ub}\sim 5$ for $\theta_{\rm rk}=30^{\circ}$), and $\theta_{\rm rk}$ is
the angle of the recoil kick with respect to the plane of the circum-black
hole disk. After the recoil, the accretion onto the black hole is expected
only from the circum-black hole disk within the radius $r_{\rm ub}$. For gas
initially within $r_{\rm ub}$, its angular momentum is rearranged by the kick
and its energy is lost due to shocks within the circum-black hole disk. The
fraction of gas within $r_{\rm ub}$ accumulated within $\sim r_{\rm{b}}$ after
the recoil kick ($f_{\rm accum}$) is found to be $f_{\rm accum}\sim
20\%$–$50\%$. Here, $\sim 20\%$ is roughly derived from the enhancement of the
surface density within $r_{\rm{b}}$ in Fig. 19 of Rossi et al. (2010), while
$\sim 50\%$ is roughly derived from the statement that all bound gas falls
within $r_{\rm{b}}$. Note that in Rossi et al. (2010), a black hole is treated
as a sink particle with a size of $0.5~{}r_{\rm{b}}$, within which gas
particles are removed from the simulation (which is described in Ponce et al.
(2012)), which presumably significantly mitigates the surface density in Fig.
19 of Rossi et al. (2010). Additionally, Rossi et al. (2010) adopted a more
centrally concentrated mass profile of $m_{\rm CBD}(<r)\propto r^{p_{\rm dm}}$
with $p_{\rm dm}=1/2$, while in the gas pressure dominated region with a fixed
opacity, the profile is less centrally concentrated, with $p_{\rm dm}=7/5$
(Haiman et al., 2009). In the case of accretion disks, the factor ($f_{\rm
inc}$) by which the gas mass and the gas density within $r_{\rm{b}}$ is
enhanced by the recoil shock is roughly estimated as

$\displaystyle f_{\rm inc}$ $\displaystyle\sim\frac{f_{\rm accum}[m_{\rm
CBD}(<r_{\rm ub})-m_{\rm CBD}(<r_{\rm{b}})]}{m_{\rm CBD}(<r_{\rm{b}})}$
$\displaystyle\sim f_{\rm accum}(r_{\rm ub}/r_{\rm{b}})^{p_{\rm dm}},$ (B3)

which is up to $\sim 60(f_{\rm accum}/0.5)$ depending on $\theta_{\rm rk}$ for
$p_{\rm dm}=7/5$.

Hence, it is numerically confirmed that the accretion rate can be enhanced by
the recoil kick in following steps. First, the gas in the annulus between
$r_{\rm{b}}$ and $r_{\rm ub}$ is perturbed by the recoil kick and it
immediately experiences shocks. Due to the shock, the gas loses its angular
momentum, it migrates inward on the dynamical timescale, and then, the gas
density in the inner regions of the circum black hole disk and the accretion
rate onto the black hole are enhanced.

Without the recoil kick, the accretion mass within $t_{\rm break}$ is $m_{\rm
CBD}(<r_{\rm acc})$, where $r_{\rm acc}$ is the radius at which gas accretes
onto the black hole on the breakout timescale, $t_{\rm break}$ ($t_{\rm
vis}(r_{\rm acc})=t_{\rm break}$, where $t_{\rm vis}(r)$ is the viscous
timescale of the circum-black hole disk at the distance $r$ from the black
hole), and is estimated as

$\displaystyle r_{\rm acc}=5\times 10^{12}~{}{\rm cm}\left(\frac{t_{\rm
break}}{0.1~{}{\rm yr}}\right)^{5/7}\left(\frac{m}{150~{}{\rm
M}_{\odot}}\right)^{1/7}$ $\displaystyle\left(\frac{{\dot{m}}c^{2}/{l}_{\rm
Edd}}{8\times 10^{4}}\right)^{2/7}\left(\frac{\alpha_{\rm
CBD}}{0.1}\right)^{4/7}\left(\frac{\kappa_{\rm CBD}}{0.4~{}{\rm
cm^{2}/g}}\right)^{1/7}$ (B4)

for the standard disk in the gas pressure dominated regime (Haiman et al.,
2009), which is expected at $r=r_{\rm acc}$ for the fiducial model but for
$f_{\rm acc}=1$ (inner regions are estimated to follow a slim disk model
(Abramowicz et al., 1988), while the disk at $r=r_{\rm acc}$ follows the
standard disk model), where $l_{\rm Edd}$ is the Eddington luminosity of the
stellar-mass black hole, $\kappa_{\rm CBD}$ is the opacity of the circum-black
hole disk, $\alpha_{\rm CBD}$ is the alpha parameter for the standard disk,
and ${\dot{m}}c^{2}/{l}_{\rm Edd}=8\times 10^{4}$ is the value in the fiducial
model with $f_{\rm acc}=1$ at $R=2.5~{}{\rm{pc}}$. After the recoil kick, the
accretion mass within $t_{\rm break}$ is

$\displaystyle\sim f_{\rm inc}m_{\rm CBD}(<r_{\rm{b}})\frac{t_{\rm
break}}{{\rm{max}}[t_{\rm break},t_{\rm vis,sh}(r_{\rm{b}})]},$ (B5)

where $t_{\rm vis,sh}(r)$ is the viscosity of the circum-black hole disk at
$r$ after the recoil shock occur in the circum-black hole disk, which can be
shorter than $t_{\rm vis}(r)$ due to the shock heating and the possible
transition to the slim disk (Abramowicz et al., 1988).

When the breakout timescale ($t_{\rm break}$) is shorter than $t_{\rm
vis,sh}(r_{\rm{b}})$, which is satisfied with $v_{\rm rk}\lesssim
300~{}{\rm{km/s}}$ for $m=150~{}{\rm M}_{\odot}$ (Eqs.B and B), the recoil
kick does not eject gas which can accrete onto the black hole within the
breakout timescale. In this case, using Eq. (B5), the accretion rate and
accordingly the averaged jet luminosity can be enhanced due to the recoil
shock by a factor of

$\displaystyle\sim f_{\rm inc}\frac{m_{\rm CBD}(<r_{\rm{b}})t_{\rm vis}(r_{\rm
acc})}{m_{\rm CBD}(<r_{\rm acc})t_{\rm vis,sh}(r_{\rm{b}})}$
$\displaystyle\sim f_{\rm inc}\frac{t_{\rm vis}(r_{\rm{b}})}{t_{\rm
vis,sh}(r_{\rm{b}})}\gtrsim f_{\rm inc}$ (B6)

where we use the rough relation of $m_{\rm CBD}(<r_{\rm{b}})/t_{\rm
vis}(r_{\rm{b}})\sim m_{\rm CBD}(<r_{\rm acc})/t_{\rm vis}(r_{\rm acc})$. In
addition, the enhancement of accretion demands that the breakout timescale be
longer than the dynamical timescale within which the shocked gas at
$r=r_{\rm{b}}$ can accrete onto the black hole, $t_{\rm
dyn,CBD}(r_{\rm{b}})<t_{\rm break}$, where

$\displaystyle t_{\rm dyn,CBD}(r_{\rm{b}})$
$\displaystyle=\left(\frac{Gm}{r_{b}^{3}}\right)^{-1/2}$ $\displaystyle\sim
8\times 10^{4}~{}{\rm{s}}\left(\frac{m}{150~{}{\rm
M}_{\odot}}\right)^{-1/2}\left(\frac{r_{\rm{b}}}{5\times 10^{12}~{}{\rm
cm}}\right)^{3/2}\,.$ (B7)

The condition ($t_{\rm dyn,CBD}(r_{\rm{b}})<t_{\rm break}$) is not satisfied
for short $t_{\rm break}$ such as in the case of GW150914-GBM, where we assume
$f_{\rm acc}\sim 1$, while that is satisfied in the case of ZTF19abanrhr with
$v_{\rm rk}\gtrsim 100~{}{\rm{km/s}}$ for $m=150~{}{\rm M}_{\odot}$.

On the other hand, when $t_{\rm vis,sh}(r_{\rm{b}})<t_{\rm break}$, satisfied
in $v_{\rm rk}\gtrsim 300~{}{\rm{km/s}}$ for $m=150~{}{\rm M}_{\odot}$, all
the mass within $r_{\rm{b}}$ accretes, and the accretion rate is enhanced by a
factor of $f_{\rm inc}m_{\rm CBD}(<r_{\rm{b}})/m_{\rm CBD}(<r_{\rm acc})$. In
this case, the accretion rate may be even reduced for $r_{\rm acc}\gg
r_{\rm{b}}$ due to the ejection of gas beyond $r_{\rm{b}}$, which is the case
when $v_{\rm rk}\gtrsim 2000~{}{\rm{km/s}}$ for $m=150~{}{\rm M}_{\odot}$ (see
Eqs. B and B). Since its effect is uncertain, we assume that the influence on
the accretion rate by the recoil shock is effectively taken into account by
enhancing or reducing $f_{\rm acc}$.

Here, the recoil velocity for GW190521 is preferred to be $\sim
300~{}{\rm{km/s}}$ or around it ($\sim 100$–$1000~{}{\rm{km/s}}$) (Abbott et
al., 2020b). Since $t_{\rm vis}(r_{\rm{b}})\sim 0.1~{}{\rm yr}$ is comparable
to $t_{\rm delay}\sim t_{\rm break}\sim 18~{}{\rm day}$ inferred for
ZTF19abanrhr, we expect that $f_{\rm acc}\sim f_{\rm inc}\gtrsim 1$, which is
predicted to be realized unless $v_{\rm rk}\gtrsim 2000~{}{\rm{km/s}}$. Since
$f_{\rm acc}$ is an uncertain quantity, when we discuss the association of
ZTF19abanrhr we treat it as a parameter whose value is adjusted to match the
observational data for ZTF19abanrhr.

Note that shocked gas, highly pressurized due to the high jet power (high
$f_{\rm acc}$), does not significantly eject the circum-black hole disk and
reduce the accretion rate onto the black hole, since the dependence of the
truncation radius of the circum-black hole disk on the accretion rate is weak
(with a power-law index of 8/53, Tagawa et al. 2022b), where the truncation
radius is the radius beyond which the circum-black hole disk is ejected by the
cocoon feedback. Also, a hollow cavity, predicted to exist around the majority
of the black holes (Kimura et al., 2021; Tagawa et al., 2022b), further
weakens this effect.

##  Appendix C Telescope Specifications

The properties of telescopes, which would be useful for detecting the
electromagnetic counterparts predicted by our model, are listed in Table 3.

Table 3The name and properties of telescopes appropriate for detecting electromagnetic counterparts  Telescope name | Photon energy | Sensitivity $[{\rm erg/s/cm^{2}}]$ | Field of view$[{\rm sr}]$
---|---|---|---
JWST |  $\sim 0.04$–$2~{}[{\rm{eV}}]$ |  $\sim 10^{-17}$–$10^{-15}$ for $t_{\rm int}\sim 10^{4}~{}{\rm{s}}$ | $\sim 10^{-6}$
WISE |  $\sim 0.05$–$0.4~{}[{\rm{eV}}]$ |  $\sim 10^{-13}$–$10^{-12}$ for $t_{\rm int}\sim 10^{4}~{}{\rm{s}}$ | $\sim 10^{-4}$
Roman space telescope (Spergel et al., 2015) |  $\sim 0.6$–$1~{}[{\rm{eV}}]$ |  $\sim 4\times 10^{-16}$ for $t_{\rm int}\sim 10^{4}~{}{\rm{s}}$ | $\sim 0.6$
ZTF (Bellm et al., 2018) |  $\sim 1.4$–$3.1~{}[{\rm{eV}}]$ |  $\sim 10^{-13}$ for $t_{\rm int}\sim 30~{}{\rm{s}}$ | 0.01
Vera Rubin (Ivezić et al., 2019) |  $\sim 1.2$–$3.9~{}[{\rm{eV}}]$ |  $\sim 10^{-15}$ for $t_{\rm int}\sim 40~{}{\rm{s}}$ | 0.003
Subaru/HSC (Aihara et al., 2018) |  $\sim 1$–$3~{}[{\rm{eV}}]$ |  $\sim 10^{-16}$–$10^{-15}$ for $t_{\rm int}\sim 10^{3}~{}{\rm{s}}$ | 0.0005
Tomo-e Gozen (Sako et al., 2018) |  $\sim 1.7$–$3.4~{}[{\rm{eV}}]$ |  $\sim 2\times 10^{-13}$ for $t_{\rm int}\sim 100~{}{\rm{s}}$ | 0.006
Chandra |  $\sim 0.2$–$10~{}[{\rm keV}]$ |  $\sim 10^{-14}$ for $t_{\rm int}\sim 2\times 10^{4}~{}{\rm{s}}$ | $6\times 10^{-5}$
XMM-Newton (Jansen et al., 2001) |  $\sim 0.4$–$3~{}[{\rm keV}]$ |  $\sim 10^{-14}$ for $t_{\rm int}\sim 10^{4}~{}{\rm{s}}$ | $8\times 10^{-5}$
HiZ-GUNDAM (Yonetoku et al., 2020) |  $\sim 0.4$–$4~{}[{\rm keV}]$ |  $\sim 10^{-11}$ for $t_{\rm int}\sim 10^{4}~{}{\rm{s}}$ | 1.2
Einstein Probe (Yuan et al., 2015) |  $\sim 0.5$–$4~{}[{\rm keV}]$ |  $\sim 3\times 10^{-11}$ for $t_{\rm int}\sim 10^{3}~{}{\rm{s}}$ | 1.0
MAXI (Matsuoka et al., 2009) |  $\sim 2$–$30~{}[{\rm keV}]$ |  $\sim 7\times 10^{-11}$ for $t_{\rm int}\sim 6\times 10^{5}~{}{\rm{s}}$ | 0.07
NuSTAR (Harrison et al., 2013) |  $\sim 10$–$30~{}[{\rm keV}]$ |  $\sim 10^{-14}$ for $t_{\rm int}\sim 10^{6}~{}{\rm{s}}$ | $3\times 10^{-5}$
FORCE (Mori et al., 2016) |  $\sim 1$–$80~{}[{\rm keV}]$ | $\sim 10^{-14}(t_{\rm int}/10^{5}~{}{\rm{s}})^{-1}$ | $10^{-5}$
Swift BAT (Barthelmy et al., 2005; Tueller et al., 2010) |  $\sim 15$–$150~{}[{\rm keV}]$ | $\sim 10^{-8}(t_{\rm int}/1~{}{\rm{s}})^{-1/2}$ | 1.4
Fermi GBM (Meegan et al., 2009) |  $\sim 8$–$4000~{}[{\rm keV}]$ |  $\sim 10^{-8}$–$10^{-6}$ for $t_{\rm int}\sim 1~{}{\rm{s}}$ | $\sim 4\pi$
INTEGRAL SPI-ACS (Winkler et al., 2003) |  $\sim 75$–$2000~{}[{\rm keV}]$ |  $\sim 10^{-7}$–$10^{-6}$ for $t_{\rm int}\sim 1~{}{\rm{s}}$ | $\sim 4\pi$

## References

  * Abbott et al. (2016a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, Phys. Rev. Lett., 116, 061102, doi: 10.1103/PhysRevLett.116.061102
  * Abbott et al. (2016b) —. 2016b, Phys. Rev. D, 93, 122003, doi: 10.1103/PhysRevD.93.122003
  * Abbott et al. (2020a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020a, Phys. Rev. Lett., 125, 101102, doi: 10.1103/PhysRevLett.125.101102
  * Abbott et al. (2020b) —. 2020b, ApJ, 900, L13, doi: 10.3847/2041-8213/aba493
  * Abramowicz et al. (1988) Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ, 332, 646, doi: 10.1086/166683
  * 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
  * Aihara et al. (2018) Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ, 70, S4, doi: 10.1093/pasj/psx066
  * 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
  * Artymowicz et al. (1993) Artymowicz, P., Lin, D. N. C., & Wampler, E. J. 1993, ApJ, 409, 592, doi: 10.1086/172690
  * Ashton et al. (2021) Ashton, G., Ackley, K., Hernandez, I. M., & Piotrzkowski, B. 2021, Classical and Quantum Gravity, 38, 235004, doi: 10.1088/1361-6382/ac33bb
  * 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
  * Barthelmy et al. (2005) Barthelmy, S. D., Barbier, L. M., Cummings, J. R., et al. 2005, Space Sci. Rev., 120, 143, doi: 10.1007/s11214-005-5096-3
  * 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
  * Barvainis (1987) Barvainis, R. 1987, ApJ, 320, 537, doi: 10.1086/165571
  * Bellm et al. (2018) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2018, Publications of the Astronomical Society of the Pacific, 131, 018002, doi: 10.1088/1538-3873/aaecbe
  * Bellovary et al. (2016) Bellovary, J. M., Mac Low, M.-M., McKernan, B., & Ford, K. E. S. 2016, ApJ, 819, L17, doi: 10.3847/2041-8205/819/2/L17
  * Berger (2014) Berger, E. 2014, ARA&A, 52, 43, doi: 10.1146/annurev-astro-081913-035926
  * Blandford & McKee (1976) Blandford, R. D., & McKee, C. F. 1976, Physics of Fluids, 19, 1130, doi: 10.1063/1.861619
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433, doi: 10.1093/mnras/179.3.433
  * Blumenthal & Gould (1970) Blumenthal, G. R., & Gould, R. J. 1970, Reviews of Modern Physics, 42, 237, doi: 10.1103/RevModPhys.42.237
  * Bromberg et al. (2011) Bromberg, O., Nakar, E., Piran, T., & Sari, R. 2011, ApJ, 740, 100, doi: 10.1088/0004-637X/740/2/100
  * Budnik et al. (2010) Budnik, R., Katz, B., Sagiv, A., & Waxman, E. 2010, ApJ, 725, 63, doi: 10.1088/0004-637X/725/1/63
  * Buonanno et al. (2008) Buonanno, A., Kidder, L. E., & Lehner, L. 2008, Phys. Rev. D, 77, 026004, doi: 10.1103/PhysRevD.77.026004
  * 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
  * Cao (2011) Cao, X. 2011, ApJ, 737, 94, doi: 10.1088/0004-637X/737/2/94
  * Cheng & Wang (1999) Cheng, K. S., & Wang, J.-M. 1999, ApJ, 521, 502, doi: 10.1086/307572
  * Connaughton et al. (2016) Connaughton, V., Burns, E., Goldstein, A., et al. 2016, ApJ, 826, L6, doi: 10.3847/2041-8205/826/1/L6
  * Connaughton et al. (2018) —. 2018, ApJ, 853, L9, doi: 10.3847/2041-8213/aaa4f2
  * Duffell et al. (2014) Duffell, P. C., Haiman, Z., MacFadyen, A. I., D’Orazio, D. J., & Farris, B. D. 2014, ApJ, 792, L10, doi: 10.1088/2041-8205/792/1/L10
  * Duras et al. (2020) Duras, F., Bongiorno, A., Ricci, F., et al. 2020, A&A, 636, A73, doi: 10.1051/0004-6361/201936817
  * Fan & Wu (2022) Fan, X., & Wu, Q. 2022, arXiv e-prints, arXiv:2212.06363.  https://arxiv.org/abs/2212.06363
  * Fan & Piran (2008) Fan, Y.-Z., & Piran, T. 2008, Frontiers of Physics in China, 3, 306, doi: 10.1007/s11467-008-0033-z
  * Farris et al. (2015) Farris, B. D., Duffell, P., MacFadyen, A. I., & Haiman, Z. 2015, MNRAS, 447, L80, doi: 10.1093/mnrasl/slu184
  * Fouka & Ouichaoui (2011) Fouka, M., & Ouichaoui, S. 2011, ApJ, 743, 89, doi: 10.1088/0004-637X/743/1/89
  * Frail et al. (2005) Frail, D. A., Soderberg, A. M., Kulkarni, S. R., et al. 2005, ApJ, 619, 994, doi: 10.1086/426680
  * Gao et al. (2013) Gao, H., Lei, W.-H., Wu, X.-F., & Zhang, B. 2013, MNRAS, 435, 2520, doi: 10.1093/mnras/stt1461
  * Gayathri et al. (2021a) Gayathri, V., Yang, Y., Tagawa, H., Haiman, Z., & Bartos, I. 2021a, ApJ, 920, L42, doi: 10.3847/2041-8213/ac2cc1
  * Gayathri et al. (2021b) —. 2021b, arXiv e-prints, arXiv:2104.10253.  https://arxiv.org/abs/2104.10253
  * Gayathri et al. (2022) Gayathri, V., Healy, J., Lange, J., et al. 2022, Nature Astronomy, 6, 344, doi: 10.1038/s41550-021-01568-w
  * Ghisellini et al. (1988) Ghisellini, G., Guilbert, P. W., & Svensson, R. 1988, ApJ, 334, L5, doi: 10.1086/185300
  * 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. (2022) Graham, M. J., McKernan, B., Ford, K. E. S., et al. 2022, arXiv e-prints, arXiv:2209.13004.  https://arxiv.org/abs/2209.13004
  * Greiner et al. (2016) Greiner, J., Burgess, J. M., Savchenko, V., & Yu, H. F. 2016, ApJ, 827, L38, doi: 10.3847/2041-8205/827/2/L38
  * Grishin et al. (2021) Grishin, E., Bobrick, A., Hirai, R., Mandel, I., & Perets, H. B. 2021, MNRAS, 507, 156, doi: 10.1093/mnras/stab1957
  * 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
  * Hankla et al. (2020) Hankla, A. M., Jiang, Y.-F., & Armitage, P. J. 2020, ApJ, 902, 50, doi: 10.3847/1538-4357/abb4df
  * 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
  * Ho (2008) Ho, L. C. 2008, ARA&A, 46, 475, doi: 10.1146/annurev.astro.45.051806.110546
  * Ito et al. (2020a) Ito, H., Levinson, A., & Nakar, E. 2020a, MNRAS, 499, 4961, doi: 10.1093/mnras/staa3125
  * Ito et al. (2020b) —. 2020b, MNRAS, 499, 4961, doi: 10.1093/mnras/staa3125
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Janiuk et al. (2017) Janiuk, A., Bejger, M., Charzyński, S., & Sukova, P. 2017, New Astronomy, 51, 7, doi: 10.1016/j.newast.2016.08.002
  * Jansen et al. (2001) Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1, doi: 10.1051/0004-6361:20000036
  * Johnson & McKee (1971) Johnson, M. H., & McKee, C. F. 1971, Phys. Rev. D, 3, 858, doi: 10.1103/PhysRevD.3.858
  * Kashiyama et al. (2013) Kashiyama, K., Murase, K., Horiuchi, S., Gao, S., & Mészáros, P. 2013, ApJ, 769, L6, doi: 10.1088/2041-8205/769/1/L6
  * Kimura et al. (2021) Kimura, S. S., Murase, K., & Bartos, I. 2021, ApJ, 916, 111, doi: 10.3847/1538-4357/ac0535
  * Kitaki et al. (2021) Kitaki, T., Mineshige, S., Ohsuga, K., & Kawashima, T. 2021, PASJ, 73, 450, doi: 10.1093/pasj/psab011
  * Kozłowski (2016) Kozłowski, S. 2016, ApJ, 826, 118, doi: 10.3847/0004-637X/826/2/118
  * Levin & Beloborodov (2003) Levin, Y., & Beloborodov, A. M. 2003, ApJ, 590, L33, doi: 10.1086/376675
  * Levinson & Nakar (2020) Levinson, A., & Nakar, E. 2020, Phys. Rep., 866, 1, doi: 10.1016/j.physrep.2020.04.003
  * LIGO Scientific Collaboration et al. (2015) LIGO Scientific Collaboration, Aasi, J., Abbott, B. P., et al. 2015, Classical and Quantum Gravity, 32, 074001, doi: 10.1088/0264-9381/32/7/074001
  * Liska et al. (2020) Liska, M., Tchekhovskoy, A., & Quataert, E. 2020, MNRAS, 494, 3656, doi: 10.1093/mnras/staa955
  * Martin et al. (2019) Martin, R. G., Nixon, C. J., Pringle, J. E., & Livio, M. 2019, New Astronomy, 70, 7, doi: 10.1016/j.newast.2019.01.001
  * Matsuoka et al. (2009) Matsuoka, M., Kawasaki, K., Ueno, S., et al. 2009, PASJ, 61, 999, doi: 10.1093/pasj/61.5.999
  * Maughan & Reiprich (2019) Maughan, B. J., & Reiprich, T. H. 2019, The Open Journal of Astrophysics, 2, 9, doi: 10.21105/astro.1811.05786
  * McKernan et al. (2020) McKernan, B., Ford, K. E. S., O’Shaugnessy, R., & Wysocki, D. 2020, MNRAS, 494, 1203, doi: 10.1093/mnras/staa740
  * McKernan et al. (2019) McKernan, B., Ford, K. E. S., Bartos, I., et al. 2019, ApJ, 884, L50, doi: 10.3847/2041-8213/ab4886
  * Meegan et al. (2009) Meegan, C., Lichti, G., Bhat, P. N., et al. 2009, ApJ, 702, 791, doi: 10.1088/0004-637X/702/1/791
  * 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
  * Mori et al. (2016) Mori, K., Tsuru, T. G., Nakazawa, K., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, T. Takahashi, & M. Bautz, 99051O, doi: 10.1117/12.2231262
  * Nakar & Sari (2010) Nakar, E., & Sari, R. 2010, ApJ, 725, 904, doi: 10.1088/0004-637X/725/1/904
  * Nakar & Sari (2012) —. 2012, ApJ, 747, 88, doi: 10.1088/0004-637X/747/2/88
  * Narayan et al. (2022) Narayan, R., Chael, A., Chatterjee, K., Ricarte, A., & Curd, B. 2022, MNRAS, 511, 3795, doi: 10.1093/mnras/stac285
  * Olano (2009) Olano, C. A. 2009, A&A, 506, 1215, doi: 10.1051/0004-6361/200912602
  * 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
  * Pan & Yang (2021) Pan, Z., & Yang, H. 2021, Phys. Rev. D, 103, 103018, doi: 10.1103/PhysRevD.103.103018
  * 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. (2016) Perna, R., Lazzati, D., & Giacomazzo, B. 2016, ApJ, 821, L18, doi: 10.3847/2041-8205/821/1/L18
  * Perna et al. (2021b) Perna, R., Tagawa, H., Haiman, Z., & Bartos, I. 2021b, ApJ, 915, 10, doi: 10.3847/1538-4357/abfdb4
  * Ponce et al. (2012) Ponce, M., Faber, J. A., & Lombardi, J. C. 2012, ApJ, 745, 71, doi: 10.1088/0004-637X/745/1/71
  * Ricci et al. (2018) Ricci, C., Ho, L. C., Fabian, A. C., et al. 2018, MNRAS, 480, 1819, doi: 10.1093/mnras/sty1879
  * 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. (2020) Romero-Shaw, I., Lasky, P. D., Thrane, E., & Calderón Bustillo, J. 2020, ApJ, 903, L5, doi: 10.3847/2041-8213/abbe26
  * 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
  * Rybicki & Lightman (1979) Rybicki, G. B., & Lightman, A. P. 1979, Radiative processes in astrophysics (Wiley and Sons, Incorporated, John)
  * Sako et al. (2018) Sako, S., Ohsawa, R., Takahashi, H., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, ed. C. J. Evans, L. Simard, & H. Takami, 107020J, doi: 10.1117/12.2310049
  * Samsing et al. (2022) Samsing, J., Bartos, I., D’Orazio, D. J., et al. 2022, Nature, 603, 237, doi: 10.1038/s41586-021-04333-1
  * Sapir et al. (2011) Sapir, N., Katz, B., & Waxman, E. 2011, ApJ, 742, 36, doi: 10.1088/0004-637X/742/1/36
  * Sapir et al. (2013) —. 2013, ApJ, 774, 79, doi: 10.1088/0004-637X/774/1/79
  * Savchenko et al. (2016) Savchenko, V., Ferrigno, C., Mereghetti, S., et al. 2016, ApJ, 820, L36, doi: 10.3847/2041-8205/820/2/L36
  * Secunda et al. (2019) Secunda, A., Bellovary, J., Mac Low, M.-M., et al. 2019, ApJ, 878, 85, doi: 10.3847/1538-4357/ab20ca
  * Sironi et al. (2013) Sironi, L., Spitkovsky, A., & Arons, J. 2013, The Astrophysical Journal, 771, 54, doi: 10.1088/0004-637x/771/1/54
  * Soldi et al. (2014) Soldi, S., Beckmann, V., Baumgartner, W. H., et al. 2014, A&A, 563, A57, doi: 10.1051/0004-6361/201322653
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757.  https://arxiv.org/abs/1503.03757
  * Spitkovsky (2008) Spitkovsky, A. 2008, ApJ, 673, L39, doi: 10.1086/527374
  * Stalevski et al. (2017) Stalevski, M., Asmus, D., & Tristram, K. R. W. 2017, MNRAS, 472, 3854, doi: 10.1093/mnras/stx2227
  * Stone et al. (2017) Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, MNRAS, 464, 946, doi: 10.1093/mnras/stw2260
  * 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. (2022a) Tagawa, H., Kimura, S. S., Haiman, Z., Perna, R., & Bartos, I. 2022a, in prep