Shock Cooling and Breakout Emission for Optical Flares Associated with Gravitational-wave Events

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

Due to the gas-rich merger environment, a key signature of the AGN channel is the possibility of electromagnetic emission accompanying the GW signal from the merger (Bartos et al. 2017;Stone et al. 2017).To explore this possibility, electromagnetic follow-up observations have been carried out for many of the mergers, with nine counterpart candidates suggested so far, including seven optical flares (Graham et al. 2020(Graham et al. , 2023) ) and two gamma-ray flares (Connaughton et al. 2016;Bagoly et al. 2016).
Recently, several studies have investigated the electromagnetic emission from a variety of transients emerging from AGN disks.Many studies (Zhu et al. 2021c,a;Perna et al. 2021a;Yuan et al. 2021;Wang et al. 2022;Lazzati et al. 2022;Ray et al. 2023) focused on the radiation from gamma-ray bursts, while Perna et al. (2021b) and Zhu et al. (2021b) discussed the electromagnetic signatures expected from accretion-induced collapse of neutron stars and white dwarfs, respectively.Yang et al. (2021), Xin et al. (2023), and Prasad et al. (2023) studied the properties of tidal disruption of stars by BHs, while Grishin et al. (2021) investigated supernova explosions, and Bartos et al. (2017) and Stone et al. (2017) estimated the electromagnetic emission produced by thermal radiation and/or outflow from circum-binary disks in AGN disks.
Several recent studies have also investigated whether transients from BHs mergingin AGN disks could explain the optical flare, ZTF19abanrhr (Graham et al. 2020), associated with the BH merger GW190521 detected in GWs.McKernan et al. (2019) discussed emission from shocks caused by collision between gas bound to the merger remnant and unbound gas after recoil kicks due to anisotropic radiation of GWs.Graham et al. (2020) assessed the net luminosity and timescales for gas accretion induced by recoil kicks.de Mink & King (2017) considered flares emerging from shocks in a circum-BH disk due to recoil kicks.Kimura et al. (2021), Wang et al. (2021a,b), and Chen et al. (2023a), respectively, considered thermal and non-thermal emission from bubbles and bubble evolution around BHs formed by strong outflows considering continuous and episodic super-Eddington accretion.Wang et al. (2021a) further considered emission from shocks emerging due to interactions of Blandford-Znajek jets (Blandford & Znajek 1977) launched from accreting BHs to the broad line regions, Rodríguez-Ramírez et al. (2023) considered free-free and bound-free emission from gas shocked due to interaction of the jets and AGN disk gas, and Tagawa et al. (2023a) estimated gamma-ray, neutrino, and cosmic-ray emission from internal shocks in the jets.Ashton et al. (2020), Palmese et al. (2021), Calderón Bustillo et al. (2021), andMorton et al. (2023) estimated the association significance of ZTF19abanrhr to GW190521.
In this paper, we develop an emission model based on the scenario proposed by Tagawa et al. (2022Tagawa et al. ( , 2023b) ) and discuss whether or not emission based on this scenario can explain some of the optical flares reported in Graham et al. (2023).Tagawa et al. (2022) indicated that a Blandford-Znajek jet can be produced from BHs embedded in AGN disks and investigated its influence on the AGN disk structure.Due to the high pressure of the shocks emerging around the jet, a cavity is created around the BHs.Just before the jet breaks out of the AGN disk, photons in the shocked gas begin to escape.These photons can be observed as breakout emission (e.g.Nakar & Sari 2010), whose properties have been investigated in Tagawa et al. (2023c).The BHs can maintain the jets even after they break out from the AGN disk, as long as there is leftover circum-BH disk gas.Once this is depleted, the BHs can no longer power the jets.This is then followed by an inactive phase which lasts until gas is replenished onto the BH, and the jet is launched again, with the cycle hence repeating.
In the case that BHs merge in the cavity while accreting (upper panel of Fig. 1), Tagawa et al. (2023b) (hereafter Paper I) predicted that electromagnetic emission is often produced associated with the BH merger.This is because the jet direction can be reoriented following a merger ( § 2.1), and strong shocks and emission are produced soon after the jet reorientation.Paper I investigated the properties of breakout emission from a jet head associated with BH mergers, and found that this model can explain various properties, including the luminosity, delay time, duration, and color of the electromagnetic transients ZTF19abanrhr, GW150914-GBM, and LVT151012, as well as why the transients began brightening only after the merger.Zhu (2023) and Zhou & Wang (2023) also estimated the neutrino emission from the breakout of the jets produced associated with BH mergers.In this paper, we additionally consider the shock cooling emission, which is produced in a subsequent adiabatic expansion phase of the shocked gas (e.g.Arnett 1980;Morag et al. 2023 for supernovae and Nakar & Piran 2017;Gottlieb et al. 2018;Hamidani & Ioka 2023 for gamma-ray bursts), in addition to the breakout emission, which is produced in an early phase of the shocked gas (see Fig. 1 for a schematic picture).We present the properties of the emission in both phases, and we then apply this model to the flares reported in Graham et al. (2023).We further discuss how to test the proposed emission model in future observations.The rest of this paper is organized as follows.In § 2, we describe a model for producing electromagnetic flares associated with GW emission and a way to constrain physical properties of the flares from the observations using this model.We present our results in § 3, discuss how to test the model in § 4, and summarize our conclusions in § 5.

METHOD
First we describe the model itself.We then specialize to discuss how to derive the model parameters from the observed properties of the flares, that is the delay time (t delay ), the duration (t duration ), the luminosity of the flare (L obs ), the merger remnant mass (m BH ), the SMBH mass (M SMBH ), and the AGN luminosity (L AGN ).In our analysis, we assume that the shocks produced by collisions between the jet and the AGN gas are characterized by non-relativistic regimes.This is because in both the breakout and the shock cooling emission, flares with the delay time and the duration of ≳ 10 day (Table 1) are usually characterized by this regime (Paper I, A schematic picture of breakout and cooling envelope emission from shocks produced by the interaction of AGN gas and a jet launched from a BH accreting gas in an AGN disk.Since the jet is reoriented at the time of the merger of two BHs ( § 2.1), electromagnetic emission is also produced after the merger.
Tables 2 and 3).As possible processes for explaining the properties of the optical flares, we consider the breakout emission from the jet head (the breakout emission scenario) and shock cooling emission from the cocoon (the shock cooling emission scenario).Note that we use the Shakura-Sunyaev model (Shakura & Sunyaev 1973) for the accretion disk in the shock cooling emission scenario, and the Thompson disk model (Thompson et al. 2005) in the breakout emission scenario.This is because the position of the BH from the central supermassive BH (SMBH) is sub-parsec for the former (Table 2) and a few parsec for the latter case (Table 3), and the mechanisms of angular momentum transfer and the disk properties are likely different in the two regions, with the Thompson disk model better suited than the Shakura-Sunyaev one to describe the outer disk.

Shock formation
Shocked gas is responsible for both the breakout and the cooling emission, and hence we begin by describing the process of shock formation.
We first describe the accretion rate onto BHs in an AGN disk, which is evaluated by the Bondi-Hoyle-Lyttleton rate.For a BH embedded in a cool AGN disk, the Bondi-Hoyle-Lyttleton radius (r BHL ) is large, and usually exceeds the scale height of the AGN disk (H AGN ) and the Hill radius (r Hill ).Accounting for the geometrical limitation of the capture regions by the shear motion and the vertical height of the AGN disk, the capture rate of gas by the BH is given by (e.g.Tanigawa & Tanaka 2016;Choksi et al. 2023), where ρ AGN is the gas density and c s,AGN is the sound speed of the AGN disk at the position of the BH (R = R BH , where R is the distance from the SMBH), v BH is the velocity of the BH with respect to the local motion of the AGN disk, v sh = r w Ω is the shear velocity at the capture radius r w = min(r BHL , r Hill ), Ω = (GM SMBH /R 3 BH ) 1/2 is the angular velocity of the BH, r h = min(r w , H AGN ) is the capture height, G is the gravitational constant, and f c ∼ 10 is a normalization constant (Tanigawa & Watanabe 2002).In the second line of Eq. (1), we assume v BH < c s,AGN and v BH < v sh .Note that r w (c 2 s,AGN + v 2 sh ) 1/2 ≈ r Hill H AGN Ω is used to derive the right hand side, which is approximately satisfied regardless of whether c s,AGN is larger or smaller than v sh .By considering the reduction or enhancement with respect to the Bondi-Hoyle-Lyttleton rate, we parameterize the fraction of the accretion rate onto the BH ( ṁ) over the Bondi-Hoyle-Lyttleton rate ( ṁBHL ) by f acc = ṁ/ ṁBHL as in Paper I. For example, low f acc may be predicted due to winds from an accretion disk with a super-Eddington rate.On the other hand, recent simulations suggest that the conversion to wind is moderate (Kitaki et al. 2021) for accretion flows in which the circularization radius (where gas is circularized after being captured by a BH) is much larger than the trapping radius (within which photons are advected to a BH without escaping), as is the case for BHs embedded in an AGN disk.In addition, the accretion rate onto a BH in a cavity during the active phases is estimated to be lower by a factor of a few compared to that without a cavity (Tagawa et al. 2022).
From rapidly accreting and spinning BHs in an AGN disk, a Blandford-Znajek jet is expected to be launched, as outlined in Appendix A.1 of Tagawa et al. (2022).The jet kinetic luminosity (L j ) is proportional to the mass accretion rate onto the BH ( ṁ) as where c is the speed of light, and η j is the jet conversion efficiency, which is approximated by η j ∼ a 2 BH for a magnetically dominated jet (Tchekhovskoy et al. 2010;Narayan et al. 2021), a BH is the dimensionless spin of the BH, and a BH ∼ 0.7 for the merger remnants (The LIGO Scientific Collaboration et al. 2021).
At a BH merger (t = 0, where t is the time from the merger), the jet direction is reoriented and can collide with unshocked AGN gas in the following ways.Once two BHs merge, the BH spin direction is reoriented if the angular momentum direction of the merging binary is misaligned with respect to the spin directions of the merging BHs.This is expected for mergers in an AGN disk due to frequent binary-single interactions (Tagawa et al. 2020a) and/or inhomogeneity of AGN disks.Since the jet is injected in the direction of the BH spin, if the jet is not aligned with the circum-BH disk due to a strong jet power (Polko & McKinney 2017), the jet propagates in the direction of the BH spin, and can collide with AGN gas.Even if the jet aligns with the angular momentum direction of the circum-BH disk due to magnetic interactions on average, the jet can precess by interacting with magnetic fields while the angular momentum directions of the BH spin and the circum-BH disk are still misaligned with one other (Liska et al. 2018(Liska et al. , 2021)).Due to the precession, the jet can collide with unshocked gas after merger during the first precession cycle (after that, the opening angle of the cavity becomes wider than that of the precession).The other possibility is that once BHs merge, a merger remnant BH receives a recoil kick in an almost random direction due to anisotropic radiation of GWs.Then shocks form in a circum-BH disk within ≲ 10 13 cm (Tagawa et al. 2023b), and shocked gas accretes onto the remnant with the angular momentum direction being modified as a result of shocks (Rossi et al. 2010).Due to magnetic interactions (Liska et al. 2021), the jet is then aligned with the angular momentum direction of the circum-BH disk, which is in turn misaligned with respect to the jet direction before the merger, and the jet can therefore collide with AGN gas.
Once the jet collides with unshocked AGN gas around a BH, two shocks form: a forward shock propagating in the AGN disk and a reverse shock in the jet.The region sandwiched by the two shocks is called the jet head.In the jet collimated regime considered in our work, the dimensionless velocity of shocked gas in the jet head at the shock breakout is estimated as (Bromberg et al. 2011) where θ 0 is the opening angle of the injected jet, and t break is the delay time between the production of the jet and the shock breakout, and is roughly given by t break ∼ 3 5 Here β FS ∼ (7/6)β h is the dimensionless velocity of the forward shock of the jet head, and f corr is the correction factor for the delay time due to the inclination of the jet and the geometry of the cavity.Ignoring geometrical corrections due to the existence of the cavity, we assume f corr = min[1/cosi, 1/θ 0 sini], considering both the cases in which the shocked gas (cocoon) breaking out of the AGN disk is from its head and from its sides, where i is the angle between the jet and the angular momentum direction of the AGN disk.With this prescription, f corr ranges between ∼ 1-1/θ 0 .
After the shock breakout, radiation is produced.Early emission is characterized by the breakout emission as described in § 2.3 and Paper I, while later emission is characterized by the shock cooling radiation described in § 2.2.Since breakout emission is associated with the shock propagation, both non-thermal and thermal emissions are expected.On the other hand, we only consider thermal emission from shock cooling, and neglect any non-thermal emission from expanding ejecta.The latter can also be bright if additional strong shocks form when the ejecta collide with the interstellar medium.We use the non-thermal component of the breakout emission and the thermal shock cooling emission to model the flares reported by Graham et al. (2023), since these are bright in optical bands.Note that the peak frequency of the thermal breakout emission from the jet head falls above the X-ray bands, which is constrained by the duration of the flare (e.g.Eqs.19 and 25 below), and hence it cannot reproduce the optical flares.

Physical model
In the shock cooling emission, photons diffuse and are released from deep inside the shocked gas, as observed in supernovae (Arnett 1980).To present the properties of this emission, we first describe the evolution of the shocked gas.
At the breakout of the shocked gas (cocoon) at t = t break , the thermal energy of the cocoon is roughly given by where β c ≃ β h θ 0 is the dimensionless velocity of the cocoon (Bromberg et al. 2011), and m BO is the mass of the cocoon at the breakout and is roughly given by considering the cylindrical shape of the cocoon with the height of H AGN f corr and the radius of H AGN f corr θ 0 .After the shock passage (even before the shock breakout), the internal energy of the shocked gas is converted to kinetic energy due to the expansion caused by the radiation pressure.Afterwards, the shocked ejecta expands nearly spherically with velocity v ej ∼ β c c.As the ejecta expands to size R, the optical depth of the spherically expanding ejecta declines as (e.g.Sapir & Waxman 2017; Nakar & Piran 2017) where κ ej is the ejecta's opacity.We adopt the Thomson scattering opacity of κ ej ∼ 0.4 g/cm 2 assuming ionized gas.If the ejecta is initially optically thick, as assumed in our fiducial model, photons deep inside the ejecta can be diffused out once the optical depth is reduced to at the time when the corresponding radius is R diff = t diff v ej and the density Due to the adiabatic expansion, the thermal energy at the diffusion radius is reduced by a factor of ∼ (R 3 BO /R 3 diff ) γ−1 where γ = 4/3 is the adiabatic index for the radiation pressure dominated gas, R BO = (V BO /4π) 1/3 is the typical size, and V BO is the volume of the cocoon at the shock breakout.Hence, the luminosity at τ ej = c/v ej is related to e BO as Using R BO in Eq. ( 11), the unperturbed AGN density at the position of the BH (R = R BH ) can be estimated via Furthermore, from the AGN density, the inflow rate of the AGN disk at R = R BH can be calculated via assuming an alpha viscosity with parameter α, and where Ω is the orbital angular velocity of the BH around the SMBH.Here, the AGN luminosity (L AGN ) is related to the inflow rate as where ṀSMBH is the accretion rate onto the SMBH, η rad is the radiation efficiency, and f cons ≡ ṀSMBH / Ṁinflow ≤ 1 is the fraction of the inflow rate at R = R BH over the accretion rate onto the SMBH.The radiation temperature is determined as follows.At τ ej = c/v ej , the energy density of the radiation within the ejecta shell is The blackbody temperature of the radiation is then given by where a is the radiation constant.Note that the radiation pressure dominates the gas pressure at t ≲ t diff for the models considered in this paper.Since the flares in Graham et al. (2023) were found by the Zwicky Transient Facility (ZTF), and the observed optical luminosity (L obs ) can be roughly estimated as the total observed energy of the flare in the optical band divided by the duration (provided in Table 3 of Graham et al. 2023), we pose that L obs needs to match the luminosity in the ZTF bands as where ν up = c/(400 nm) and ν down = c/(700 nm) are the upper and lower limits for the frequencies observed by the r and g bands of the ZTF.Note that the luminosities in the r and g bands of the ZTF are not directly derived from Eq. ( 17).To calculate them (panels c and d of Fig. 3), we adjust the limits of the integral in this equation to cover the frequency range of the relevant band.
2.2.2.Derivation of model parameters Using Eqs. ( 1)-( 17), we can then determine the 17 and Ṁinflow given the input parameters, θ 0 , f acc , a BH , α, η rad , f corr , and f cons , and the observed properties, L obs L AGN , and t diff .For example, Combining Eqs. ( 1)-( 14), the velocity of the cocoon is calculated as where f jet/BHL ≡ f c f acc η j is a parameter related to the jet power, and A and δ are variables depending on the velocity of the cocoon.For β c /θ 0 < 1, A = (35/18) 6 θ −3 0 and δ = 1/2, and otherwise (as long as the jet is in the collimated regime) A = θ 0 and δ = 1/8.To determine β c using Eq. ( 18), L SC needs to be derived using Eqs.( 15)-(17).To consistently solve these equations, we calculate β c in an iterative way using Newton's method.Then, by incorporating β c into Eqs.(1)-( 17), the other parameters can also be determined.
We further adjust f jet/BHL so that the scale height of the AGN disk is equal to the height expected for the Shakura-Sunyaev disk, given α.Such disk structure is expected to be realized in regimes where the Toomre parameter (Eq.20) becomes greater than 1.

Breakout emission 2.3.1. Physical model
Long before photons deep inside the ejecta escape (t ∼ t break + t diff ), at around the time that the shock arrives to the surface of the ejecta (t ∼ t break ), bright breakout emission is released.Since breakout emission of the jet head is associated with the shock propagation, both non-thermal and thermal emission are expected.Such non-thermal emission can explain the optical flares found by Graham et al. (2020), as investigated in Paper I. In the following we describe how to reconstruct the model properties for the breakout emission.
Photons inside the shock start to diffuse out from the AGN disk and the breakout emission begins to be released when the photon diffusion time from the shock front, t diff,sh ∼ d 2 edge κ sh ρ AGN /c, becomes equal to the dynamical timescale of the shock, t dyn,edge ∼ d edge /β FS c, where d edge is the thickness of the AGN disk above the shock, and κ sh is opacity of the shocked gas.By equating these timescales, the thickness at the breakout is given by d edge,BO ∼ 1/(β FS κ sh ρ AGN ), and the duration of the emission from a breakout shell is where we adopt κ sh ∼ 0.4 g/cm 2 considering the ionization of gas by photons released from the shocks (Nakar & Sari 2010).
If we assume that the AGN disk is gravitationally unstable at the position of the BH, which is expected for the delay time and the duration of the flares (Paper I), the density of the AGN disk is related to the position of the BH through (e.g.Thompson et al. 2005) where Q is the Toomre parameter.Note that AGN disks at the BH positions are predicted to be Toomre unstable (Table 2, e.g.Haiman et al. 2009).
In the case of breakout emission, R BH is found to be large ( § 3.2).At such large scales, efficient transfer of angular momentum of the AGN gas is required for SMBH accretion (Thompson et al. 2005).Following Thompson et al. (2005), we assume that the inflow rate of the AGN disk is parameterized by 2.4.Evolution In this section, we describe the evolution of the luminosity and the temperature for thermal emission from the cocoon in the breakout and shock cooling emission phases.Referring to previous studies (Arnett 1980;Nakar & Sari 2010;Piro et al. 2021;Morag et al. 2023), we assume that the luminosity evolves as where t ′ = t − t break , t pl = 1/(β 2 c cκ ej ρ AGN ) is the duration of emission from a breakout shell in the cocoon, t sph ≃ R BO /β c c is the transition between the planar and spherical geometries of the breakout shell, t dyn = H AGN /β c c is the dynamical timescale, and n is the power law slope of the vertical AGN gas density profile at the height at which photons begin to break out.The second and third rows correspond to the luminosity in the planar phase (before the shocked gas doubles its radius by expansion) for the breakout emission, the fourth row corresponds to that in a spherical phase (after the shocked gas doubles and before the shock cooling emission phase), and the fifth row corresponds to that in a shock cooling emission phase.
Similarly following previous studies, we assume that the temperature evolves as where is the black body temperature of the breakout shell, n T is the power-law index for the temperature evolution in the spherical phase, and n T = −(18.48n 2 + 20.69n + 6)/[(1.19n+ 1)(22.32n+ 17)] for an expanding spherical ejecta (Nakar & Sari 2010).The validity of these analytical formulae has been tested by recent numerical simulations (e.g.Piro et al. 2021;Morag et al. 2023).Since Eqs. ( 23) and ( 24) are formulae for an expanding ejecta with a spherically symmetric evolution, some modifications are required in the cases of emission from an expanding ejecta with a cylindrical shape like a cocoon, especially at around the transition between the planar

Table 1
The properties of the ZTF flares suggested to be associated with the GW events.In each pair, the flare is located within the 90% credible volume of and detected within the delay time of 200 days from the GW event.Here, M SMBH is the SMBH mass, m BH is the merged remnant mass, t duration is the duration of the flare, t delay is the delay time between the detection of the GWs to and the flare, Lop is the optical luminosity of the flare, and L AGN is the Eddington ratio of the host AGN.For J053408.41 + 085450.6, the SMBH mass is not constrained and a fiducial mass of 10 8 M ⊙ is used as in Graham et al. (2023).
Pair LIGO/Virgo alart ID Name of AGN log 10 and spherical phases.For simplicity, we determine n and n T so that the luminosity and the temperature smoothly evolve at t ′ = t diff , respectively.

Parameters
In this section we describe the fiducial values for the model parameters, and the observed properties of the flares used in the modellings.
As fiducial values, we set the opening angle of the injected jet to θ 0 = 0.2 (e.g.Hada et al. 2013Hada et al. , 2018;;Berger 2014), the radiation efficiency to η rad = 0.1, the correction factor for the delay time to f corr = 3 (which roughly corresponds to the median value considering the cavity with the aspect ratio of ∼ 1 and isotropic jet directions), the consumption fraction of the inflow rate to f cons = 1, the alpha-viscous parameter to α = 0.1, and the angular momentum transfer parameter in the outer regions of the AGN disk to m = 0.1.For the computation of the shock cooling emission we adjust the value of f jet/BHL ≡ f c f acc η j so that the scale height of the AGN disk is equal to the height expected for the Shakura-Sunyaev disk ( § 2.2), while for the breakout emission we set it to f jet/BHL = 10 assuming f c = 10 (Tanigawa & Watanabe 2002), η j = 0.5, and f acc = 2 considering moderate enhancement of accretion due to shocks caused by recoil kicks (Appendix B, Tagawa et al. 2023b).
When modeling the flares with breakout emission, we assume that the flare duration (t duration ) corresponds to the exponential decay time (t e ) in Graham et al. (2023).To derive the delay time (t delay ), we identify the day at which the flare luminosity peaks using digitizer1 .Due to the difficulty of identifying the peak, t delay contains uncertainties.We adopt M SMBH and m BH from Tables 3  and 4 of Graham et al. (2023).The optical luminosity is derived by means of the total observed energy in the optical band divided by the sum of the rise time (t g ) and the decay time (t e ) presented in Graham et al. (2023).We calculate the luminosity of the AGNs (L AGN ) by inferring the flux of the AGNs at ∼ 4000 Å from Fig. 6 of Graham et al. (2023), and estimating the luminosity distance from the redshift of the GW events assuming a value for the Hubble constant of 67.8 km/s s −1 Mpc −1 , for the matter density today of 0.24, and for the cosmological constant today of 0.74 (Planck Collaboration et al. 2016), and adopting the bolometric correction factor of 5 (Duras et al. 2020).The values for the observed quantities adopted in this paper are listed in Table 1.Note that the delay time (t delay ) corresponds to t break , and the duration of a flare (t duration ) corresponds to t diff,BO .
Conversely, when modeling the flares with shock cooling emission, we assume that the duration of a flare (t duration ) corresponds to t diff .The delay time between a GW event and an optical flare (t delay ) is also on the order of ∼ t diff given t break ≪ t diff , while we do not use the delay time (t delay ) to constrain the model parameters.Note that t delay ∼ t duration expected in this scenario is roughly consistent with the properties of the observed flares (table 1).We assume that the observed optical luminosity of the flare (L op ) corresponds to the luminosity in the ZTF bands (L obs ).
We note that in our model, the physical properties are uniquely determined.This is because we fixed several input parameters as detailed above, and we do not directly use the observational data points for parameter fitting.If instead we allowed the input parameters to vary, we would introduce degeneracies among several input parameters.The variations of the values of the fixed input parameters can affect physical properties of the model especially in the breakout-emission scenario, while less affect them in the shock cooling emission scenario (see § 4.1 and § 4.2).We believe that such simple prescriptions are useful to understand typical properties expected from the model and consider possible tests of the model below.

Shock cooling emission
Table 1 shows the list of possible pairs of associations between the GW events and the electromagnetic flares reported in Graham et al. (2023), along with their observed properties.Seven electromagnetic flares and twelve pairs for the possible associations are reported.The number of the pairs ( 12) is larger than that of the flares (7), Table 2 The inferred model parameters for the shock cooling emission scenario.R BH is the distance from the SMBH to the BH, βc is the expansion velocity of the cocoon, h AGN ≡ H AGN /R BH and ρ AGN are the aspect ratio and the density of the AGN disk at the position of the BH, respectively, T BB and λ peak are the radiation temperature and wavelength of the shock cooling emission, respectively, L jet is the jet power, f jet/BHL is a parameter related to the jet power, t break is the breakout timescale of the jet, and t diff,CBO is the diffusion timescale for the breakout emission from the cocoon.Table 3 The inferred model parameters for the breakout emission scenario.Ṁinflow / ṀEdd is the Eddington ratio for the gas inflow at the position of the BH, f op/kin is the fraction of the optical luminosity of the flare over the kinetic power of the jet, f BH/NSC is the fraction of R BH over the size of the nuclear star cluster, and h AGN /h AGN,TQM is the ratio of the aspect ratio of the AGN disk derived from our model to that derived assuming the model in Thompson et al. (2005).In the last column, the value of the correction factor for the delay time, at which because some flares can be associated with more than one GW events.Hence, at least five pairs are false associations.Note that in this paper we do not analyze the two gamma-ray flares possibly associated with GW events due to their significantly different properties (Connaughton et al. 2016;Bagoly et al. 2016), while they can be explained by the thermal breakout emission from the jet head in relativistic regimes (Tagawa et al. 2023b).
Table 2 shows the distribution of the model parameters when the properties of the observed flares (Table 1) are modeled with shock cooling emission.
The parameter f jet/BHL is widely distributed depending on the pairs.f jet/BHL ∼ 1-10 in pairs 4, 11, and 12 is roughly expected for the Bondi-Hoyle-Lyttleton accretion, while f jet/BHL ∼ 10-60 in pairs 1-3, 5, and 7-10 can be realized by the enhancement of accretion by shocks due to recoil kicks at merger (Paper I).However, f jet/BHL ∼ 3000 in pair 6 is difficult to be realized, since only up to f jet/BHL ≲ 300 is found to be feasible by considering the radial surface density profile of a circum-BH disk and appealing to the results of numerical simula-tions (see Appendix B of Paper I).On the other hand, even for pair 6, if we adopt f corr ∼ 1, f jet/BHL is reduced to ∼ 300.This is presumably because both a high accretion rate and a low inclination have similar influences on the breakout velocity, and most parameters are mainly characterized by the magnitude of the breakout velocity.Due to the variation of f corr (between ∼ 1-1/θ 0 ) and the degeneracy between f corr and f jet/BHL , the cooling model is poorly constrained or tested by the value of f jet/BHL .
The parameter R BH ranges between 0.06-0.3pc, corresponding to 900-6 × 10 4 R g , where R g = GM SMBH /c 2 is the gravitational radius of the SMBH.These locations roughly correspond to migration traps, gap forming regions, and/or slow migration regions, where BH mergers are predicted to be frequent (Bellovary et al. 2016;Tagawa et al. 2020b).
The orange lines in the left panel of Fig. 2 show the SEDs for pairs 1-12.The radiation temperature ranges between T BB = 9200-20000 K, corresponding to the peak wavelength λ peak = ch/(3k B T BB ) = 240-520 nm, where h is the Planck constant and k B is the Boltzmann con- Figure 2. The SEDs for shock cooling (orange) and preceding breakout emission (blue) from a cocoon in the shock cooling emission scenario for pairs 1-12 (left panel), and those for thermal breakout emission (brown) from a jet head in the breakout emission scenario (right panel) assuming the luminosity distance of d L = 3 Gpc.The spectral energy distributions (SEDs) for the non-thermal emission is drawn for the model adopted in Fig. 5 of Paper I as a representative example.The dotted blue, purple, red, green, dark green, and gray lines mark the sensitivity of the ZTF, HiZ-GUNDAM and Einstein Probe, Chandra and XMM-Newton, NuSTAR, ULTRASAT, and Swift BAT, respectively, as adopted in Paper I. The vertical cyan lines present the g and r bands of ZTF.Note that the shock cooling emission is too dim (≪ 10 −17 erg/s/cm 2 ) in the breakout emission scenario to be visible in the right panel.
stant.Also, the various colors observed in the optical flares (Fig. 3 of Graham et al. 2023) can be reproduced depending on the variation of the radiation temperature with changing some input parameters, such as the direction of the jet (i).For example, if we set i = 0 • (the jet being perpendicular to the AGN-disk plane) by using derived values for ρ AGN , R BH , and H AGN and the observed parameters (m BH and M SMBH ), the radiation temperature of flares is enhanced by a factor of ∼ 1.7.This is because for i = 0 • the shocked mass becomes low, and then photons can escape from an earlier phase when the radiation temperature is higher.Thus, the variety of colors can be reproduced if there are temperature variations in the thermal emission.For synchrotron emission (whose possible contribution to the flares is discussed in § 3.2), the color is related to the power-law slope of the injected electrons accelerated by the first-order Fermi process (p) as νL ν ∝ ν (−p+2)/2 , and p is presumably distributed in a narrow range for the same phenomena (Panaitescu & Kumar 2001).On the other hand, various values for the color can be realized if the breaks of the power laws in the SEDs coincidentally fall at around the ZTF bands (see the black line in the right panel of Fig. 2).Note that the distribution of the colors is currently uncertain due to the shift in the baseline of the flux and the contribution from the background AGN emission in Fig. 3 of Graham et al. (2023).If the color is actually distributed in a wide range, the flares are then easier to model by thermal shock cooling emission rather than by non-thermal breakout emission.In addition, the model can be also discriminated by constraining the spectral shape of the emission.This can be done by simultaneously observing the flares by the ZTF and the Ultraviolet Transient As-tronomy Satellite (ULTRASAT, Shvartzvald et al. 2023, Fig. 2) in the future.We also predict breakout emission from the cocoon preceding the shock cooling emission (this is the emission from the cocoon, which is different from the emission from the jet head discussed in § 3.2).The duration of the breakout emission from the cocoon, is distributed in the range t diff,CBO ∼ 10-10 5 s (the right most column of Table 2), and the temperature (Eq.25) is in a range 2×10 5 -2×10 6 K (blue lines in the right panel of Fig. 2).Although the duration of the breakout emission is shorter than than that of the shock cooling emission, it can be detected by future X-ray surveys, such as HiZ-GUNDAM (Yonetoku et al. 2020) and/or the Einstein Probe (Yuan & Narayan 2014).Fig. 3 a and b show the evolution of the luminosity and the temperature of emission for pairs 1, 5, and 11.We chose pair 1 as a fiducial example, pair 5 as an event with small m BH , and pair 11 as the one with long t duration .The lines are drawn until the phase at which τ ej = 1 is satisfied, since our model for thermal emission is not valid when τ ej < 1.
For the smaller m BH case (dashed lines in Fig. 3), the timescale for the breakout emission of the cocoon (t pl ) is longer.This is because L j is lower due to the lower Bondi-Hoyle-Lyttleton rate (Eqs. 1, 2), and β c tends to be lower for lower L j (Eqs.3, 18).To reproduce t diff , ρ AGN needs to be lower for lower β c (Eq. 9).Due to the low β c and ρ AGN (Table 2), the breakout timescale of the cocoon (t pl ∝ β −2 c ρ −1 AGN ) becomes longer.In such ), apparent magnitude in the g (mg, blue) and r (mr, red) bands, AB magnitude at a near-ultraviolet band (260 nm, black, panel c), and the color (mg − mr, panel d) for thermal emission in the shock cooling emission model for pairs 1 (solid), 5 (dashed), and 11 (dotted).We present pairs 1, 5, and 11 as a fiducial, a smaller m BH , and a longer t duration case, respectively, among which evolution timescales are derived to be different.Lines are drawn until τ ej = 1 after which our model is invalid.
case, since photons diffuse out from the shocked gas earlier, adiabatic expansion is inefficient, leading to higher temperature.
For the longer t duration model (dotted lines in Fig. 3), the evolution of the luminosity is slower in the shock cooling regime (t ′ ≳ t diff ).

Breakout emission
Table 3 shows the parameters of the breakout emission model reproducing the properties of the observed flares.In this model, the parameters are distributed in narrow ranges.The mergers are predicted to occur at R BH ∼ 1-10 pc, the shock velocity ranges within β h ∼ 0.3-0.4,the disk aspect ratio at R BH ranges within h AGN ≡ H AGN /R BH ∼ 0.001-0.02,the fraction of the optical luminosity to the kinematic power of the shock is around L op /L j ∼ 0.02-1, the fraction of the distance from the SMBH to the size of the nuclear star cluster (R NSC ) ranges within R BH /R NSC ∼ 0.2 − 0.7 assuming the empirical relation of R NSC = 8.3 pc(M SMBH /3.1 × 10 8 M ⊙ ) 0.154 (Scott & Graham 2013;Kormendy & Ho 2013;Georgiev et al. 2016), and last the inflow rate of an AGN disk in units of the Eddington rate varies within Ṁinflow / ṀEdd ∼ 0.004 − 1.2, where ṀEdd ≡ L Edd /c 2 η rad with radiative efficiency of η rad = 0.1, and L Edd is the Eddington luminosity of the SMBH.We find that the optical luminosity of the flare is lower than the kinetic power of the jet (f op/kin < 1) for all the pairs (Table 3).Considering the bolometric correction of the breakout emission in the optical band (∼ 10, Paper I), our scenario demands f op/kin ≲ 0.1.For the pairs that do not satisfy this condition, we need to use f jet/BHL > 10.Note that a higher value (f jet/BHL ≲ 100) is possible by considering the enhancement of the accretion rate due to shocks caused by recoil kicks at merger (Appendix.B of Paper I).
In the right panel of Fig. 2, the SEDs for thermal emission from the breakout of the jet head for Pairs 1-12 are displayed.As they are bright at ∼ 10 17 -10 19 Hz and the duration is long (≳ 10 6 s), they can be observed by X-ray telescopes, such as the Swift X-ray telescope (XRT, Burrows et al. 2005), Chandra, XMM Newton (Jansen et al. 2001), and the Nuclear Spectroscopic Telescope Array (NuSTAR, Harrison et al. 2013).Note that the shock cooling emission associated with the breakout emission in all the pairs is so dim (≪ 10 40 erg/s) that it is buried within the AGN variability.
The derived values of R BH /R NSC ∼ 0.2−0.7 ∼ 1−8 pc are much larger than in the shock-cooling case.The different locations are driven by the phases of the emission.In our model, the duration of flares corresponds to the timescale at which the diffusion timescale is equal to the dynamical timescale (e.g.Eqs. 10 and 19 for the shock cooling and breakout emission cases, respectively).This timescale is inversely proportional to the density of the ejecta or of the AGN disk, as well as to the square of the velocity of the ejecta or of the shocks in the shock cooling and breakout emission scenarios, respectively.As a result, since the velocities are typically comparable (Tables 2 and 3), the density of the AGN disk in the breakout emission scenario must be similar to the density of the ejecta in the shock cooling emission scenario, in order to reproduce the duration of the flares.Additionally, in the shock cooling emission scenario the density of the ejecta is much lower than the local AGN density due to adiabatic expansion (Eq.12).Thus, the AGN density is typically much lower in the breakout emission scenario, and therefore the distance from the SMBH needs to be larger.These larger distance are consistent with the scenario that BHs in nuclear star clusters are captured by AGN disks and merge with one other.Here, the migration timescale of objects in AGN disks around massive SMBHs as inferred for the flaring AGNs (Table 1) is so long that migration before BH mergers is less expected.
To check whether the values derived for the aspect ratio (h AGN ) are plausible, we calculate the aspect ratio of the AGN disk (h AGN,TQM ) adopting the model in Thompson et al. (2005), using R BH , M SMBH , Ṁinflow , and m.We find that the values of h AGN,TQM are roughly comparable to h AGN within a factor of ∼ 2 (9th column of Table 3).Such moderate differences could arise due to the varia-tion of f corr reflecting the variation of the inclination of jets.In Table 3, we list f corr at which h AGN,TQM = h AGN .This shows that h AGN,TQM = h AGN is satisfied by reasonable values for f corr around ∼ 3-5.
By comparing L AGN /L Edd in Table 1 and Ṁinflow / ṀEdd in Table 3, we can determine the value of m at which L AGN /L Edd = Ṁinflow / ṀEdd as Ṁinflow ∝ m.For the pairs 11 and 12, m ≥ 3 is required to satisfy L AGN /L Edd ≤ Ṁinflow / ṀEdd .In the other pairs (1-10), m can be as low as ≥ 0.003-0.5 to satisfy this condition.We presume that m ≲ 1 is acceptable (e.g.Hopkins & Quataert 2010;Begelman & Armitage 2023;Begelman & Silk 2023), although there are no studies constraining the possible ranges of m as far as we know.Note that the inflow rate at pc scales can be much higher than the accretion rate onto the central SMBH (e.g.Izumi et al. 2023), since a large fraction of gas is possibly consumed by star formation and outflows, as predicted in Blandford & Begelman (1999), Thompson et al. (2005), and Collin & Zahn (2008).Then Ṁinflow / ṀEdd is allowed to be much larger than L AGN /L Edd , and the required value of m is enhanced.Hence, pairs 11 and 12 are not well explained by breakout emission due to the required high values of m.On the other hand, it is not clear whether AGN disks are usually in steady state, which might complicate the comparison between the inflow rate and the AGN luminosity.

TESTS OF THE MODEL
In the following we discuss how our model can be tested by examining the distribution expected for the observable properties.From § 3.1 and 3.2 we see that values for each model parameter are distributed in a narrow range (Tables 2 and 3).We discuss whether this is because flares originating from merging BHs tend to have well-defined characteristic properties, or because of observational selection effects in the way the search was conducted by Graham et al. (2023).Flares were searched with specific ranges of parameters, and in particular with the rise time within 5-100 day, the decay time within 10-200 day, and the delay time from the GW detection in the interval 0-200 day.
In order for electromagnetic flares to be found in association with BH mergers, BHs need to typically merge in bright AGNs.This is because most BH mergers reported by LIGO/Virgo/KAGRA are found to merge at luminosity distances of several Gpc (Abbott et al. 2020).At such large distances, AGNs are easily missed unless they are bright.Indeed, the host SMBH masses for the flares reported by Graham et al. (2023) are ≳ 10 8 M ⊙ , hence so massive that AGNs are rarely missed in AGN searches at the distances of the GW events.Also, assuming a luminosity distance of d L ∼ 3 Gpc, SMBH mass of 10 8 M ⊙ , Eddington ratio of ∼ 1, bolometric correction of 5, and the fraction of the variable luminosity compared to the average luminosity ∼ 0.1, the flux of variable flares in the AGN is ∼ 2 × 10 −13 erg/s/cm 2 , which falls just above the sensitivity of ZTF ∼ 10 −13 erg/s/cm 2 .Conversely, it is difficult to find flares associated with BH mergers occurring in AGNs around less massive SMBHs through surveys for AGN variability.
We next estimate the detection rate of electromagnetic flares associated with BH mergers based on our scenario.If BH mergers actually produce electromagnetic flares as found in Graham et al. (2023), the rate of such flares is comparable to or less than the rate of BH mergers (∼ O(10) Gpc −3 yr −1 ).If flares within ≲ 3 Gpc are observable by the ZTF and all merging BHs produce electromagnetic flares, up to N BBH,3Gpc ∼ 300 flares associated with BH mergers can be detected per year.On the other hand, we have estimated that a small fraction of BH mergers (f EM/BBH ≲ 0.02) could accompany detectable electromagnetic flares in our model (see § 4.1 in Paper I for discussions).Then, the seven flares discovered by (Graham et al. 2023) is roughly consistent with the number of flares expected during LIGO/Virgo/KAGRA O3, which is estimated as N EM/GW,O3 ∼ N BBH,3Gpc f EM/BBH t O3 ≲ 10, where t O3 ∼ 1.1 is the duration of the O3 operation in the unit of year.Constraints on the frequency of flares, whose properties can be explained by emission from BHs, are hence useful to test our model.Here, note that dimmer more frequent flares would be contributed by solitary BHs (Tagawa et al. 2023c).

Shock cooling emission
From Table 2, we note that there are moderate variations in R BH and T BB in the shock cooling emission scenario.To consider a possible test of this model, we discuss its dependence on the observed properties, and the range of these properties for which the shock cooling emission scenario is inconsistent.If there are events with high or low T BB , extremely large or small R BH , or R BO /R diff ≥ 1, such flares are inconsistent with the shock cooling emission scenario.However, the dependence of the variables T BB , R BH , R BO /R diff on the ratio of the observed parameters (L AGN , L obs , m BH , M SMBH , and t diff ) over the adjustment parameters (f corr and f jet/BHL ) is similar to the dependence of H AGN /R BH on the same ratio.For example, for β c /θ 0 < 1, H AGN /R BH depends on the observed properties as An interesting test of the model is the correlation between the delay time and the duration of the flare.This is because the delay time is comparable to the duration of a flare as long as t diff ≫ t break , which is satisfied in pairs 1-12 (Table 2).To constrain the correlation coefficient e.g. with uncertainty of ≲ 0.3 by 95 percentile, more than ≳ 50 events are needed to be observed.Hence, more events will be very helpful as further diagnostics.Note that half of the pairs reported by Graham et al. (2023) are unreal as a flare is associated with a few GW events, and hence, we need to derive the correlation excluding the influence from false associations.Also, we derived the delay time using digitizer.In addition, the time dependence of the luminosity assumed by Graham et al. (2023) is different from that expected in the shock cooling emission.Thus, both the delay time and the duration suffer from significant uncertainties.If these timescales can be well constrained, the model can be tested for each event by comparing the delay time and the duration.
Another possible test is the detection of the cocoon breakout emission preceding the shock cooling emission.In pairs 1-12, the duration and temperature of the cocoon breakout emission range in the intervals ∼ 10-10 5 s and 0.4-5 keV; thus wide-field X-ray surveys, such as Einstein Probe (Yuan et al. 2015) and HiZ-GUNDAM (Yonetoku et al. 2020) will be useful to detect the early breakout emission (left panel of Fig. 2).Assuming that the luminosity of the breakout emission is similar in magnitude to the jet power, the detectable distance is estimated to be where F sen is the sensitivity of the facility.Here, the closest distance among the flares discovered by Graham et al. ( 2023) is ∼ Gpc.Since the duration of the cocoon breakout emission is ∼ 10-10 5 s in pairs 1-12 and the sensitivity of Einstein Probe (Yuan et al. 2015) and HiZ-GUNDAM is F sen ∼ 10 −11 erg/s/cm 2 for t int ∼ 10 4 s, events at the luminosity distance of ∼ Gpc can be detected if L jet ≳ 10 46 erg/s.If both the breakout emission and the shock cooling emission are detected from the same AGN, it can be a robust test of this model.
It is notable that in the shock cooling model, the color keeps evolving to redder in all pairs (Fig. 3 d).The m g − m r evolves by ∼ 0.4 in 50 days.Such mild evolution of the color can be a strong test of this model.On the other hand, for non-thermal emission from the breakout of the jet head, the color is expected to be unchanged.Note that the colors are presented for pair 3 in Fig. 2 of Graham et al. (2020), revealing almost no evolution.However, due to the contamination from the host AGN emission, any color evolution is likely diffi-cult to constrain.To derive the color evolution more precisely, observations at additional frequencies, e.g. by ULTRASAT (black lines of Fig. 3 c), would also be useful.Hence, when available, the evolution of the color will be an additional diagnostic in future observations.

Breakout emission
Next we discuss ways to test the breakout emission scenario.To do this, we first present the dependence of physical quantities on observables as For the breakout emission from the jet head, if the shocked gas becomes relativistic (β h ≳ 0.8), the probability of observing non-thermal emission is significantly reduced due to the relativistic beaming effects, the shift of the minimum energy, and the time dilution (as shown in Fig. 3 b of Paper I).Thus, if the shock is relativistic, the breakout emission is presumably not observed by current facilities, such as ZTF.Since β h is estimated to be ∼ 0.2-0.4 for the observed flares, if β h is higher by a factor of ∼ 3 compared to our estimates, flares cannot be explained by breakout emission.β h ≳ 0.8 is satisfied in all the pairs if events with delay time of t delay ≲ 2 day are found.Note that this is not compensated by f corr and f jet/BHL (as in the shock cooling emission scenario).This is because β h is reduced only by a factor of ∼ 1.1 by enhancing f corr to the maximum value (∼ 1/θ 0 ).Also, if f jet/BHL is reduced to lower β h , f op/kin becomes larger than 1, which violates another requirement for the model.If associations between the flares and GWs are due to random coincidence, the delay time is expected to be distributed uniformly in the range of 0-200 day.Assuming a uniform distribution, we can test the breakout emission scenario at the ∼ 1 σ and ∼ 2 σ levels after discovering ≳ 100 and ≳ 300 events, respectively, by checking if there are events with t delay ≲ 2 day.If we find optical flares with t delay ≲ 2 day, the breakout emission scenario is disfavored.
If t duration is enhanced by one order of magnitude, the jet power is reduced by a similar factor, and the fraction of the optical luminosity to the jet power exceeds one for Pairs 1-3, 6, 11, and 12, and then, the breakout emission scenario becomes inconsistent for these pairs.However, the enhancement of L j due to long t duration can be compensated by f jet and f corr up to about a factor of ∼ 10, which is limited by the requirement for β h (≲ 0.8) as discussed above.Hence, the model is not well tested by t duration .
Currently, the properties of the observed flares satisfy the conditions required for the breakout emission scenario.To be consistent with this scenario, the delay time should not be shorter than ∼ 2 day, and the color of flares should be distributed in a narrow range as discussed in § 3.1.Thus, to test whether there are flares with properties being inconsistent with the breakout emission scenario, more events will need to be observed.

CONCLUSIONS
In this paper we have presented the properties of emission from shocks emerging from collisions between AGN gas and a jet launched from a merger remnant BH in an AGN disk.Our model includes the evolution throughout the shock breakout and subsequent cooling emission phases.We then applied this model to the candidate flares reported in Graham et al. (2023).Our results are summarized as follows.
1. We fit the characteristic features of all of the events with both emission processes.Both processes could fit each observation, suggesting that such fits may themselves be insufficient to rule out the AGN origin for the flares.The reconstructed parameters might then be indicative of the selection algorithm determining the false alarm rate of associations.
2. While both processes could be made consistent with the observed events with appropriate parameter selection, we found that the implied merger distance from the central SMBH is markedly different for the two processes.Specifically, shock cooling emission can explain the observed properties if the mergers happen R BH ∼ 0.06-0.3pc from the SMBH, which is consistent with the locations of AGN-assisted merger models (e.g.Tagawa et al. 2020b).On the other hand, breakout emission would require a much larger distance of R BH ∼ 1-8 pc to explain the observed flare duration.This may be possible for mergers in AGNs with high SMBH masses, in which migration is inefficient (e.g.Perna et al. 2021b).
3. Follow-up observations could help further constrain the reconstructed parameters of the events.X-ray observations would have to be made prior to the optical detection, which likely requires future widefield surveys.In addition, follow-up observations determining the spectral evolution of the electromagnetic flares would be important.
H.T. was supported by the National Key R&D Program of China (Grant No. 2021YFC2203002) and the National Natural Science Foundation of China (Grant No. 12173071).S.S.K. was supported by Japan Society for the Promotion of Science (JSPS) KAKENHI Figure1.A schematic picture of breakout and cooling envelope emission from shocks produced by the interaction of AGN gas and a jet launched from a BH accreting gas in an AGN disk.Since the jet is reoriented at the time of the merger of two BHs ( § 2.1), electromagnetic emission is also produced after the merger.

Figure 3 .
Figure3.The evolution of the luminosity (panel a), temperature (panel b), apparent magnitude in the g (mg, blue) and r (mr, red) bands, AB magnitude at a near-ultraviolet band (260 nm, black, panel c), and the color (mg − mr, panel d) for thermal emission in the shock cooling emission model for pairs 1 (solid), 5 (dashed), and 11 (dotted).We present pairs 1, 5, and 11 as a fiducial, a smaller m BH , and a longer t duration case, respectively, among which evolution timescales are derived to be different.Lines are drawn until τ ej = 1 after which our model is invalid.
duration [day] t delay [day] Lop [erg/s] L AGN /L Edd AGN , L obs , m BH , M SMBH , and t diff ) are observed.Thus, in this model, once H AGN /R BH is adjusted to the value expected in the Shakura-Sunyaev model, T BB , R BH , and R BO /R diff are then characterized by realistic values, and the distribution of the observables becomes difficult to use as a further test of the model.
Then, if we adjust f corr and f jet/BHL so that H AGN /R BH is consistent with a Shakura-Sunyaev disk, these parameters (T BB , R BH , and R BO /R diff ) also fall in a range of possible values (similar values to those derived in Table 2), even if events with wide ranges of observed pa-rameters (L