Optical Emission Model for Binary Black Hole Merger Remnants Travelling through Discs of Active Galactic Nuclei

Active galactic nuclei (AGNs) have been proposed as plausible sites for hosting a sizable fraction of the binary black hole (BBH) mergers measured through gravitational waves (GWs) by the LIGO-Virgo-Kagra (LVK) experiment. These GWs could be accompanied by radiation feedback due to the interaction of the BBH merger remnant with the AGN disc. We present a new predicted radiation signature driven by the passage of a kicked BBH remnant throughout a thin AGN disc. We analyse the situation of a merger occurring outside the thin disc, where the merger is of second or higher generation in a merging hierarchical sequence. The coalescence produces a kicked BH remnant that eventually plunges into the disc, accretes material, and inflates jet cocoons. We consider the case of a jet cocoon propagating quasi-parallel to the disc plane and study the outflow that results when the cocoon emerges from the disc. We calculate the transient emission of the emerging cocoon using a photon diffusion model typically employed to describe the light curves of supernovae. Depending on the parameter configuration, the flare produced by the emerging cocoon could be comparable to or exceed the AGN background emission at optical, and extreme ultraviolet wavelengths. For instance, in AGNs with central engines of ∼ 5 × 10 6 M ⊙ , flares driven by BH remnants with masses of ∼ 100 M ⊙ can appear in about ∼ [10-100] days after the GW, lasting for few days.


INTRODUCTION
A natural prediction of stellar dynamics is the clustering of thousands of stellar mass black holes (sBHs) orbiting a few parsecs around the central supermassive black holes (SMBHs) of galaxies (Bahcall & Wolf 1977;Morris 1993;Miralda-Escudé & Gould 2000;Hailey et al. 2018).In active galactic nuclei (AGNs), these sBHs can receive torque forces bringing their orbits aligned to the plane of the AGN disc (McKernan et al. 2012(McKernan et al. , 2018)), and migrate into certain radial regions, denoted as migration traps (Bellovary et al. 2016;Secunda et al. 2019).Under this scenario, such migrations induce close encounters among compact objects, thus making AGNs discs plausible candidates for hosting the binary black holes (BBH) coalescences detected by the LIGO-Virgo-Kagra (LVK) experiment (The LIGO Scientific Collaboration et al. 2021;Bartos et al. 2017).
Remnants from non-symmetrical binary mergers can be born with recoil or kick velocities of ∼[100 -1000] Km s −1 (Campanelli et al. 2007a,b;Zlochower & Lousto 2015;Varma et al. 2022).In typical AGNs, these kicked remnants can be retained by the AGN potential well, having then the chance to undergo further encounters with compact objects.Thus, AGNs could also host BH hierarchical growth, i.e., BHs growing through binary mergers where one or both components are the remnants of a previous coalescence.Hierarchical ★ E-mail: juancr@cbpf.brmerging in AGNs is a viable channel to explain the observed anticorrelation among the mass ratio and effective spin of BH mergers (Santini et al. 2023) as well as mergers with components of masses ≳ 50 M ⊙ (Chatziioannou et al. 2019;Kimball et al. 2020;Abbott et al. 2020a), whose existence challenges standard stellar evolution models (Woosley 2017;Abbott et al. 2020b).Hierarchical growth is also a possible channel for the formation of intermediate-mass black holes (IMBHs) (McKernan et al. 2012(McKernan et al. , 2014)), a class of BHs with mass in the range [10 2 -10 5 ] M ⊙ , thought to be the seed of SMBHs, and showing currently poor observational evidence compared to their stellar mass and super-massive counterparts (Greene et al. 2020).
An appealing property of the AGN merger channel, not present in other BBH merger scenarios, is the ability to generate multimessenger emissions.Contrary to gravitational waves from compact binaries involving neutron stars, BBH mergers are not expected to produce electromagnetic (EM) counterparts by themselves.However, BBHs coalescing nearby high-dense media, like the thin discs of AGNs, may produce multimessenger emission (e.g., GWs and radiation) due to the interaction of the binary and/or the remnant with the gas of the disc (McKernan et al. 2019;Graham et al. 2020;Kimura et al. 2021;Wang et al. 2021;Tagawa et al. 2023b,a).
Among the GWs measured by the LIGO-Virgo experiment, the event GW190521 has been of particular interest in the context of mergers assisted by AGNs.This GW event produced a BH remnant of ∼ 142 M ⊙ , the heaviest BH measured by an LVK observation so far, from a merger with inferred components of 85 and 55 M ⊙ (Abbott et al. 2020a).Such massive binary components are unlikely of stellar collapse origin (Woosley 2017) and more naturally explained by the hierarchical merger channel.Moreover, an optical flare in the AGN J1249 + 3449 at redshift  = 0.438 was claimed EM counterpart to GW190521 by Graham et al. (2020), who interpreted the EM flare as originated in a hyper-accretion episode onto the merger remnant.This GW event has one of the largest localisation volumes on the sky among the events measured by LVK so far (Palmese et al. 2021), which makes the multi-messenger association a subject of ongoing debate (Ashton et al. 2021;Palmese et al. 2021;Chen et al. 2022;De Paolis et al. 2020).If established, the aforementioned and future EM counterparts to GWs can serve as direct probes of AGN discs, could be employed to derive cosmological constraints (e.g., Haster 2020 ), as well as an independent method for measuring the Hubble constant (Abbott et al. 2021;Gayathri et al. 2020;Bom & Palmese 2023).Further theoretical and observational analyses are then timely to solidly localise EM signals associated with GW events.
This paper presents a new predicted EM signature triggered by BBH mergers in AGNs.We consider a BBH merger occurring outside and nearby the plane of an AGN thin disc, whose kicked remnant eventually plunges and traverses the disc.We study one of the jet cocoons (Bromberg et al. 2011) driven by the remnant within the disc, focusing on the case when the cocoon propagates quasi-parallel to the disc.The jet cocoon eventually emerges outside the disc in the form of a non-relativistic outflow.We then calculate the emission produced after the emerging cocoon expands and let the thermal photons escape.The proposed emission scenario is motivated by previous works that consider the extraction of disc material as a viable mechanism to explain thermal and non-thermal flares from AGNs (Ivanov et al. 1998;Pihajoki 2016;Valtonen et al. 2019;Rodríguez-Ramírez et al. 2020).
The present EM counterpart scenario is also motivated by second or higher-generation mergers in a hierarchical sequence which could explain mergers detected by LVK with components more massive than predicted by stellar evolution models (e.g., GW190521 Abbott et al. 2020b andGW170729 Abbott et al. 2019).The components of such binary systems have then been born in previous coalescences whose recoils perturb the alignment of the remnants with the disc plane.Thus, mergers of second or higher generations could occur within or outside the disc and here we consider the latter case (different to the previous works of McKernan et al. 2019;Graham et al. 2020;Kimura et al. 2021;Wang et al. 2021;Graham et al. 2023;Tagawa et al. 2023b, that consider mergers within the disc).
The paper is organised as follows.Section 2 describes the proposed emission scenario and derives an analytic model for the emitted spectrum.In Section 3, we discuss the parametric space of interest and present the spectral energy distributions and light curves predicted by the model.The visibility and temporal signature of the flare is analysed in Section 4. Finally, we summarise and discuss our findings in Section 5.

THE EMERGING COCOON SCENARIO
Highly spinning BHs with kick velocities of ∼ [100-1000] Km s −1 can result after coalescences of non-equal mass binaries (Campanelli et al. 2007a,b;Zlochower & Lousto 2015).In addition, spinning BHs accreting magnetised gas can launch jets with 100% efficiency, or higher, due to extraction of spin energy as described by the Blandford-Znajek mechanism (Blandford & Znajek 1977;McKinney et al. 2012;Kaaz et al. 2022).Motivated by the aforementioned theoretical predictions, here we analyse the situation of a spinning BH remnant of BBH coalescence that enters the dense region of an AGN thin disc and launches relativistic jets while travelling within it.We focus on BH remnants from mergers of second or higher generations, which can be viable explanations for high-mass remnants inferred by LVC detections, like GW170729 (Abbott et al. 2019) and GW190521 (Abbott et al. 2020b).Since the components of a second-generation merger were produced in previous coalescences, these components were likely born with post-merger kicks that set them outside of the disc for a significant fraction of their orbital period, as we estimate in Appendix A. When such kicked remnants then form a new binary system, the orbit of the latter could also be misaligned with the disc plane.This motivates us to explore the case of a second (or higher) generation merger occurring outside the AGN disc.Thus, in the present analysis we consider a BBH coalescence taking place at a distance  from the AGN SMBH, at a vertical distance from the disc mid-plane  0 , being ℎ d () <  0 < , with ℎ d () the disc semi-thickness at the radial distance .A kicked remnant of mass  • is produced at coalescence with instantaneous kick velocity  k forming an angle  k with respect to the disc normal, as depicted in Figure 1A.After the coalescence, the BH remnant meets the disc surface in the (source frame) period (1) Once the remnant enters the disc, it undergoes Bondi-Hoyle-Lyttleton (BHL) accretion (Bondi & Hoyle 1944;Lee et al. 2014;Lora-Clavĳo et al. 2015) of the disc material.Such accretion drives relativistic jets which take to form a period of after the BH enters the disc.In equation (2), is the so-called BHL radius, and   is a factor that depends on the magnetisation of the external medium.In this work we consider BH remnants launching jets with the fiducial efficiency of 200%.According to general relativistic magneto-hydrodynamical (GRMHD) simulations of wind accretion performed by Kaaz et al. (2023), jets with 200% efficiency can be launched by BHs travelling within a gaseous environment of plasma beta  ∼ 10.Such efficient jets take a period of ∼ 2 HL / k to form, as can be seen in Figures 3 and 4 of Kaaz et al. (2023).Thus, we consider the fixed value of   = 2 through this paper.
As the BH trajectory proceeds, the jets propagate within the disc inflating bipolar cocoons (Bromberg et al. 2011) of shocked disc material that eventually meets the edge of the disc (see Figure 1B).Hence, we consider the jet cocoons as the channel through which the BH remnant transports mass and energy outside the disc.Here we refer to the cocoon material that emerges from the disc as the "emerging cocoon".
When the cocoon meets the disc boundary (defined by the disc semi-height), its jet-like morphology is no longer preserved due to the drastic drop of gas density beyond the disc semi-height ℎ d1 .At this point, the cocoon material prefers to expand laterally on the side where the disc density drops, as depicted in Figure 1C.Here we focus on emerging cocoons driven by jets propagating quasi-parallel to the disc plane.We note that compared to jets propagating in a quasiperpendicular direction, quasi-parallel jets have the chance to sweep up more disc material and their cocoons break out the disc edge with a slower flow velocity along the disc's vertical direction (quasiperpendicular jets would drive faster and more jetted-like outflows).Thus, we approximate the emerging cocoon as a quasi-spherical, non-relativistic expanding flow.When breaking out the disc, this outflow could produce a brief high-energy EM transient due to the acceleration of non-thermal particles at the outflow expansion front.Here we are, however, focused on the thermal emission produced by photons that diffuse within the emerging cocoon and emanates from its outer surface.
The emerging cocoon is composed by swept-up disc material collected while the BH travels within the disc.At the same time, the jet that creates the cocoon takes the period given by the equation (2) to form.Therefore, the production of the emerging cocoon is constrained by the thickness of the disc where the remnant plunges.In the present analysis, we then define as a necessary condition for the production of emerging cocoons.
In the next subsection, we calculate the mass and energy transported by the cocoon outside the disc on the observer side and in Subsection 2.2, we assess the expansion of the emerging cocoon and its thermal emission.

The jet cocoon propagation within the thin disc
We consider the energy of the emerging cocoon as equivalent to the energy injected by the travelling BH through one of its jets and within the period when the BH undergoes BHL accretion of the disc material.For analytic purposes, we assume uniform gas density within the disc (with a thickness of semi-height ℎ d ) and drastically smaller density outside.We then estimate the energy stored in the emerging cocoon as where  j is the BH jet power and Δ bh is the time interval since the jet is formed until the BH reaches the AGN disc boundary (see Figure 1).The time interval of jet energy injection then is being the first term in the RHS of equation ( 6) the time that the BH takes to cross the disc of thickness 2ℎ d , and the second term the time needed for the jet to form (see equation 2).Motivated by the results of Kaaz et al. (2022), we take the power of one of the BH jets as which corresponds to a total jet power with an efficiency of 200% (the BH launches two jets), being the accretion rate onto the travelling BH.The RHS of equation ( 8) represents a fraction  acc of the BHL accretion rate onto a massive particle travelling within a non-magnetised medium with sound speed  s and gas density  d .Through this paper we adopt the fixed value of  acc = 0.1, which is consistent with a magnetised medium of  ∼ 10, according to Kaaz et al. (2023).
We estimate the mass of the emerging cocoon as the mass of the material enclosed within the jet cocoon just before breaking out the disc.The cocoon encloses a mixture of jet-ejected material plus swept-up disc material.We assume the jet-ejected gas to be much more diluted compared to the disc material, and that the mass of the cocoon is always dominated by the swept-up mass when the cocoon meets the disc edge.Then we estimate the mass of the cocoon as follows: where we approximate the volume of the cocoon as a cylinder of height  H and radius  c , being Δ  the time interval since the jet is formed (see equation 2) until the cocoon reaches the AGN disc edge (see Figure 1B).To obtain the height and radius of the cocoon, we employ the formalism of Bromberg et al. (2011), for the case of a relativistic jet propagating through a uniform medium.Then, we note that the period Δ c is related to the disc half-thickness and the remnant's velocity as: where  H and  c are the height and the cylindrical radius of the cocoon, parameterised as: We adopt the solutions derived by Bromberg et al. (2011) to obtain the speeds of growth  H and  c (in units of the speed of light) for the height and radius of the cocoon, respectively: with  ∼ 0.7, and  0 the jet opening angle at the base.To obtain  H ,  c and hence the mass within the cocoon (equation 9), one can first solve equation ( 10) to obtain Δ  and then evaluate equations ( 11).
We adopt the standard thin disc model of Shakura & Sunyaev (1973) (hereafter SS), to parameterise the properties of the disc (such as the semi-height ℎ d , density  d , and speed of sound  s ) as a function of the distance to the SMBH , and its mass  s and accretion rate  s .The cocoon solution given by equations ( 12)-( 13) corresponds to the "collimated jet" regime (Bromberg et al. 2011).This is a suitable description for the cocoon evolution as long as L < 1.If on the other hand L >  −4/3 0 , one should adopt the mathematical solutions corresponding to the "uncollimated jet".Combining equations ( 7), (8), and ( 14) we estimate In this work we explore EM counterparts from BH remnants of  • ≲ 200 M ⊙ and  k ≳ 200 Km s −1 .In addition, we find that solutions to equation ( 10) give in general Δ c ≫ 0.1 day.Thus, the parameter configurations considered here lead in general to L ≪ 1 and we then consider the collimated jet solution only throughout this work.

The expansion of the emerging cocoon
The outflow of disc material that emerges from the disc is a mixture of matter and photons that produce an EM flare after expanding enough to let escape the thermal photons.For calculating such emission, here we model the evolution of the emerging cocoon as equivalent to an expanding sphere of mass  0 , of initial uniform density  0 , and total energy  0 , similarly to a supernova remnant.
We take the mass and initial volume of the expanding sphere as the mass and volume of the jet cocoon just before breaking out the disc (see equation 9 and related text in the previous section), and hence we take the initial density and radius as  0 =   / 0 and  0 = [3 0 /(4)] 1/3 , respectively.We assume the total energy  0 (given by equation 5) to split into kinetic and thermal energies during the outflow expansion: respectively and we assume energy equipartition setting  k = 1/2.The emerging cocoon is radiation pressure dominated when breaking out the disc.Thus the initial pressure  0 , and temperature  0 can be related as being  r the radiation constant.Simultaneously, the initial pressure and volume can be related to the thermal energy (equation 17 ) as being  a = 4/3 the adiabatic index appropriate for a radiation pressure dominated gas.Since the density of the environment outside the disc is negligible compared to the outflow density, we consider a free expansion for the outflow outer radius  =  0 +  0  ′ , with the constant radial velocity We calculate the emission of the emerging cocoon as that of a supernova (SN) remnant in its free expansion phase following (Arnett 1980(Arnett , 1996;;Chatzopoulos et al. 2012).In this approach, the emission of the spherical expanding plasma is produced by photons arriving by diffusion at the outflow surface.The peak of bolometric luminosity occurs when the diffusion timescale equals the dynamic timescale.Following the approach of (Arnett 1996;Chatzopoulos et al. 2012), the diffusion time is where  is the gas opacity taken as constant, and  = 4 3 /9.Considering the outflow dynamic time as  h = / 0 , the condition  d =  h occurs at Henceforth, the bolometric luminosity evolves as (Arnett 1996;Chatzopoulos et al. 2012): with   = √︁ 2 d,0  h,0 ,  d,0 = /( 0 ), and  h,0 =  0 / 0 .We take  as the electron opacity  T =  T /( e  u ), being  T the electron scattering cross-section,  u the atomic mass constant,  e = 2/(1 + ), and we use  = 0.85 as the hydrogen mass fraction.We then calculate the spectrum emitted by the emerging cocoon as black-body radiation of effective temperature being  SB the Stefan-Boltzmann constant.Thus, the flux density of radiation at the time  and frequency  in the observer frame (at Earth) is being  and  L the source red-shift and luminosity distance, respectively, and   ′ is the black-body spectral radiance of temperature  eff .

AGN+EMERGING COCOON EMISSION PROFILES
To investigate the wavelengths at which the outflow flare discussed in Section 2 can be observed, we compare its spectrum with the emission of the hosting AGN.To obtain a particular emission profile, we first specify assumed values for  (redshift of the source),  L (luminosity distance),  s (SMBH mass),  s (SMBH accretion rate), , (disc viscosity parameter), and  (distance from the SBMH where the merger occurs).Based on these parameters, we derive through the standard disc model of SS, the properties of the local environment where the remnant BH interacts, namely  d (gas density),  d (temperature), and ℎ d (disc half thickness).Then, specifying the parameters of the remnant, namely  • (BH mass)  k (kick velocity),  k (kick angle relative to the disc normal), and  0 (jet opening angle), we obtain the thermal emission from the emerging cocoon as described in Section 2.
Given the chosen values of  s and  s , we assess the emission of the hosting AGN as  ,bg = (  n / ref ) ,ref , where  ,ref is the average AGN spectrum profile given by the blue or cyan points in Figure 7 (Ho 2008),  ref is the reference luminosity of this spectrum at  ref = 4400 Å, and  n is a normalisation factor that depends on the mass and accretion rate onto the SMBH and that we define following Tagawa et al. (2023b) as: In this normalisation, we use the value of  c = 3, leading to a background AGN spectrum consistent with the calculated disc emission.
Given the values of  s ,  s , and , we calculate the disc emission by integrating the black body radiation of the disc surface (see e. g., Frank et al. 2002), from the innermost stable circular orbit of a non-spining SMBH, 6 g , to 10 4  g .Motivated by the location of the AGN J124942.3+344929,claimed as the host of the first EM counterpart to a BBH event (Graham et al. 2020), throughout this paper we use the fiducial values of  = 0.438 and  L = 1734.166Mpc, for the redshift and luminosity distance of the source, respectively2 .The location where the BH remnant might interact with the disc is so far not well constrained, thus we explore multiple locations within the range of  = [1000 − 8000]  g ( g =   s / 2 ), which are of the order of the migration trap location discussed in Bellovary et al. (2016).
In order to assess the disc properties at the above radii, we employ the disc model of SS and consider accretion rates within [0.01, 0.1]  Edd with the Eddington accretion rate defined as  Edd = 1.39 × 10 18 [/ ⊙ ] g s −1 .Since thin accretion discs can be gravitationally unstable at parsec scales, we evaluated the radius within which the disk is stable against self-gravity varying the SMBH mass between 10 6 and 10 10  ⊙ (Collin-Souffrin & Dumont 1990).We find that SS discs are stable for SMBHs with   ≲ 10 8 M ⊙ at the radii and accretion rates considered here.Thus, in this paper, we explore EM counterparts restricted to SMBHs ranging [5 × 10 6 − 5 × 10 7 ] M ⊙ .This choice for the SMBH mass range is also be motivated by the number density of SMBHs at the local Universe, which is higher for masses ∼ [10 6 − 10 7 ] M ⊙ than for masses ≳ 10 8 M ⊙ (Ueda et al. 2014).
We assume that recoils from previous mergers slightly perturbed the alignment with the disc of the binary components discussed here, and thus, we consider the merger to occur outside the disc.The vertical location of the merger  0 is then restricted to the maximum displacement from the disc mid-plane that the binary components attained due to the kicks of their previous coalescences.In AGNs with SMBHs of masses ≳ 5 × 10 6 M ⊙ , remnants of mergers occurring within the disc at ∼[1000 -8000]  g from the SMBH and with kick velocities ≲ 1000 Km s −1 , are retained by the AGN gravitational potential.Such kicked BHs can reach a maximum vertical displacement from the disc mid-plane ranging about ∼ [5-100] times of the disc semi-height ℎ d , as we estimate in Appendix A. This maximum −displacement can be reduced due to dynamical friction acting every time the kicked BH crosses the thin disc.Detailed calculation of the remnant's orbit considering the interacting with the AGN disc is beyond the scope of the present work.Thus, we consider the initial position of the merger considered here (a second or higher order merger generation) as a free parameter.To illustrate the solutions derived from the present multi-messenger scenario, we consider the values of  0 /ℎ d = 5 and 20, which are consistent with the maximum −displacements derived in Appendix A.
Assuming that the jet propagates quasi-perpendicular to the BH trajectory, we choose the small angle of  k = 8 • , since we are focusing on jets propagating quasi-parallel to the disc plane (see Section 2).For the jet opening angle, we use the fixed value of  0 = 15 • in all calculations in this work.We explore emission profiles corresponding to BHs with masses and kick velocities in the ranges of  M ⊙ and [100-1000] Km s −1 , respectively, motivated by potential measures of BHs by the LVK experiment (Abbott et al. 2020a) as well as by numerical models of non-symmetrical BBHs mergers (Zlochower & Lousto 2015).
In the following subsection, we discuss the time delay at the Earth frame for the appearance of the EM counterpart and in Subsections 3.2 and 3.3, we present spectral energy distributions (SEDs) and light curves (LCs), respectively, of the flare emission.

The flare starting time
In the multi-messenger scenario discussed here, a BBH coalescence produces a GW signal followed by an EM flare starting after a time delay Δ pGW .This time delay is computed as: which comprises the temporal periods in which (i) the BH remnant enters the disc Δ ent (equation 1), (ii) the jet forms within the disc Δ j (equation 2), (iii) the jet cocoon reaches the disc boundary Δ c (equation 10), and (iv) the emerging cocoon produces its maximum bolometric luminosity Δ max (see equation 22).In equation ( 29), the factor (1 + ) accounts for the time dilation due to the red shift of the source.We illustrate in Figure 2 the delay Δ pGW and its components as a function of the kick velocity.This example corresponds to a remnant of 150 M ⊙ , from a coalescence at  0 = 5ℎ d and  = 5000 g from an SMBH of 3 × 10 7 M ⊙ accreating at 0.03  Edd .We note that in the limit of  0 ≫ ℎ d the Δ ent component dominates in the RHS of equation (29).In this limit, we can approximate The other three components represent a subdominant contribution to the delay time Δ pGW and we simply approximate them as the time that the remnant spends within the disc ≈ 2ℎ d / k .Thus in the  0 ≫ ℎ d limit, we approximate Δ pGW as: In Figure 3, we display the delay Δ pGW (equation 29) as a function of the kick velocity, for different values of the mass of the remnant (indicated by the curve thickness), and different SMBH masses and accretion rates.Upper and lower panels are obtained assuming  s = 6 × 10 6 and 3 × 10 7 M ⊙ , respectively, whereas blue and red curves correspond to 0.03 and 0.1  Edd , respectively.In each panel, the upper and lower sets of curves correspond to  0 = 20 and 5ℎ Km s −1 travelling through discs around SMBHs of 5 × 10 6−7 M ⊙ drives flares starting ∼[3-300] days after the GW signal, accordingly We note that the delay Δ pGW increases for larger mass and accretion rate onto the SMBH.The delay  pGW is insensitive to the mass of the BH remnant when  0 is relatively large, as expected (see equation 30).

Spectral energy distributions
Figure 4 displays spectral energy distributions (SEDs) corresponding to the differential luminosity of the emerging cocoon, the AGN disc, and the AGN background emission as measured at the Earth.Blue curves are the outflow SEDs at different time steps after the GW event, where the first SED curve is calculated at the time when the photon diffusion and the dynamical time of the emerging cocoon coincide (see equation 22).The black dashed curves plot the stationary disc spectra, and grey curves plot the AGN emission.All curves in this figure were obtained assuming a remnant BH of 100 M ⊙ with kick velocity of 800 Km s −1 , interacting at  = 5000 g and  0 = 5ℎ d with the disc.Panels on the left and right correspond to an SMBH of 6 × 10 6 and 3 × 10 7 M ⊙ , respectively, whereas upper and lower panels were obtained assuming accretion rates of 0.03 and 0.1  Edd (onto the SMBH), respectively.The examples of SED curves of Figure 4 illustrate some general features of the EM counterpart flare described in Section 2. We note that for the relatively low accretion rate of 0.03  Edd the flare can outshine or produce emission comparable to the AGN at NIR, optical and EUV wavelengths.For the larger AGN accretion rate of 0.1  Edd , the flare emission peaks at lower frequencies.The examples of SED profiles in Figure 4 then hints that the EM counterpart studied here could be well observed in EUV and optical, and perhaps in NIR.In the following subsection, we describe the starting times and duration In each panel, the curves are calculated using a fixed distance  as indicated.Upper and lower panels are obtained assuming an SMBH of 6 × 10 6 and 3 × 10 7 M ⊙ , respectively.Blue and red curves correspond to the accretion rates of 0.03 and 0.1, in Eddington units, respectively.The sets of upper and lower curves in each panel are obtained assuming  0 = 20 and 5ℎ d , respectively.
of such outflow flares by exploring their LC profiles at NIR, optical, and EUV wavelengths for different masses and kick velocities of the BH remnant.

Light-curves
We derive light-curve (LC) profiles as a superposition of AGN plus flare emissions observed at Earth.For simplicity, we consider the AGN background emission (grey curves in Figure 4) as stationary, and we do not model the rise of the flare.Instead, we add the flare to the background emission at the time when the emerging cocoon emits its maximum bolometric luminosity (see Subsection 2.2).The obtained LC profiles are shown in Figure 5 as flux per unit frequency where we illustrate examples of LCs at NIR, optical, and EUV frequencies.The assumed GW event occurs at  pGW = 0.The curves constantly start in time (stationary AGN emission) eventually displaying a flux jump followed by a temporal evolution that converges to the initial background emission.The left and right panels in Figure 5 are obtained assuming a remnant BH of 75 and 150 M ⊙ , respectively, whereas plots in the upper and lower panels assume an SMBH of 5 × 10 6 and 10 7 M ⊙ , respectively.In each panel, LCs of fluxes at NIR, optical, and EUV frequencies are plotted in brown, green and blue colours, respectively, being the upper curve triads associated to the SMBH accretion rate of 0.1  Edd , and the lower triads to the 0.01  Edd rate.Different line styles are associated to different remnant kick velocities, as indicated.
The kick velocity has noticeable effects on the flare starting time, amplitude, and duration.For progressively larger values of  k in the range of [600 -1200] Km s −1 , we note that in general: (i) the flare starting time decreases, (ii) the flare amplitude decreases, (iii) the decaying flare period increases.We also note that in AGNs with  larger accretion rates onto the SMBH, the emerging cocoon produces flares with larger decaying times.The mass of the BH remnant has the opposite effect: BHs with larger masses produce flares with shorter decay times.

OBSERVATIONAL FEATURES
AGNs typically exhibit continuum variability on timescales from weeks to years, currently of unknown origin.The emission variations of non-blazar AGNs are observed on the order of 10% of their base emission in optical-UV wavelengths (Vanden Berk et al. 2004;Yu et al. 2022).This AGN intrinsic variability can challenge the identification of the flares derived in this work as GW counterparts.
In addition, different astrophysical processes can also produce flaring and enhancements in the AGN emission, such as variations in the accretion onto the SMBH (Ross et al. 2018), magnetic reconnection in the vicinity of the SMBH (de Gouveia Dal Pino et al. 2010;Scepi et al. 2021), supernova and kilonova events in the AGN disc (Grishin et al. 2021;Zhu et al. 2021), and tidal disruption events (Rees 1988;Chan et al. 2019).
We are then interested in studying the conditions where the flares of the emerging cocoon exceed the AGN typical variability within the shortest possible time lags, as these solutions could be better identified as GW counterparts.In the following subsections, we explore the parameter space of the present model to study the behaviour of the amplitude, time lag, and duration of the predicted flares.We define as distinguishable those flare solutions representing more than 50% of the hosting AGN, or equivalently, magnitude differences of |Δ| ≳ 0.5 mag in optical bands3 .

Fractional excess and time delay
The observed frequencies where the flare could be better distinguished over the AGN background can vary depending on the parameter configuration of the system (i.e., AGN plus the BH remnant).
Here we study the visibility of the emerging cocoon, considering the fractional excess   = (  −  ,bg )/ ,bg which measures the intensity of the flare flux   relative to the AGN background emission  ,bg .In Figure 6, we display   as a function of the mass of the remnant for different values of  s (SMBH mass),  k (kick velocity), M BH = 75 M Figure 5. LC profiles obtained from the emerging cocoon model described in Section 2. The curves are flux densities at selected emission wavelengths, namely NIR(J band, 1235 nm, brown), optical (g-band, 464 nm, green), and extreme ultraviolet (UVw2 band, 193 nm, darkcyan).The LCs were obtained assuming an SMBH of 5 × 10 6 M ⊙ (upper) and 10 7 M ⊙ (lower).As indicated, the left and right panels assume BH remnants of different masses and different curve styles correspond to solutions with different BH kick velocities.The upper and lower curves in each panel correspond to accretion rates of 0.1 and 0.01 (in Eddington units), respectively, of the AGN central engine.The AGN is assumed located at a redshift of  = 0.438.The timescale on the −axis indicates the time elapsed, in the observer frame, after the assumed GW event.
and  (distance from the SMBH).We calculate   with   at  max (see subsection 2.2), which corresponds to the maximum bolometric luminosity of the flare.We consider BH kick velocities within the range of [250 -1000] Km s −1 , and the three wavelengths chosen for the light curves of Figure 5, representing the NIR, optical, and EUV domains.The end of the curves in this figure indicates that the condition (4) for the production of an emerging cocoon is no longer satisfied, according to the analysis of Section 2. We indicate with the red horizontal lines the values of   = -0.5, 0, and 0.5, which correspond to flares emitting 50%, 100% and 150% of the background emission.
The panels on the left in Figure 6 explore the effect of different SMBH masses, considering the fixed location of  = 5000 g .We note that in AGNs with SMBHs of 5 × 10 6 M ⊙ , BH remnants should have masses lower than ∼ 100 M ⊙ to produce detectable flares, and such flares can be produced by remnants with masses as low as ∼ 30 M ⊙ .As the mass of the SMBH increases, the mass of the BH remnant should be larger to produce detectable flares.In an AGN with a SMBH of 5 × 10 7 M ⊙ , for instance (left lower panel of Figure 6), remnants should be heavier than ∼ 50M ⊙ to produce flares exceeding the background emission.The right panels of Figure 6 are obtained with the fixed SMBH mass of 10 7 M ⊙ and explore the effect of different locations  where the remnant interacts with the disc.We note a trend where for larger values of , the fractional excess   is enhanced.We also note that fluxes calculated at the EUV wavelength (blue curves) are, in general, more visible than their optical and NIR counterparts (blue and brown curves, respectively).This can be seen by comparing the fractional excess   among the curves corresponding to the same kick velocity (i.e., with the same line style).We note that the EM counterpart can be seen at NIR provided that the BBH merger occurs at distances of  ≳ 8000 g , in AGNs with SMBHs ≳ 5 × 10 7 M ⊙ .Nevertheless, at those radii the Shakura & Sunyaev disc model employed here could be no longer valid (see Section 3).
In Figure 7, we show the predicted time delay Δ pGW in the observer frame after the GW event (thick curves) for the appearance of the flare and its duration Δ dec (thin curves) as a function of the mass of the remnant using the same parameter configurations considered in Figure 6 (for  s ,  s ,  k , and ).The time delay for the appearance of the flare is calculated with equation (29).We take the flare duration as the elapsed time among the flare appearance and when the initial flux diminishes a factor of  = 2.718.The curves in the upper panels of Figure 7 explore different SMBH masses using the fixed value of  = 5000 g , whereas the lower panels explore different values of  with the fixed SMBH mass of 10 7 M ⊙ .We note that flares are faster as  s and  are shorter.For instance, detectable flares lasting as short as ∼ 2 − 5 days can be produced at radii of ∼ 2000 g , (see left-lower panel in Figure 7).On the other hand, relatively prominent flares like the ones on the left-lower panel of Figure 6, are relatively slower.This strong flares are produced by remnants interacting at  ≳ 8000 g and take ∼[100-300] days ([500-1000] days) to appear for  0 = 5ℎ d ( 0 = 20ℎ d ), having durations of ∼ [10-30] days.Moreover, we note that Δ pGW and Δ dec increase for larger  s and .Finally, we note that Δ pGW do not vary significantly with the mass of the remnant, and Δ dec decreases as the mass of the remnant is larger.

Optical light-curve and detectability
The flares produced by the emerging cocoon can be comparable to or exceed the emission of the hosting AGN at optical bands and last for a few days to a few weeks within the parameter space of interest (see Figures 6-7).Thus, the flares predicted in this work can be detected by optical time-domain surveys capable of capturing transients of a few days or longer, such as the ongoing Zwicky Transient Facility (ZTF; Bellm et al. 2019), the forthcoming Vera C. Rubin Observatory (LSST; Ivezić et al. 2019), and possible systematic searches using DECam (Bom et al. 2023;Morgan et al. 2020).The detectability of the flares also depends on the distance of the source as well as the sensitivity of the observing instrument.In Figure 8, we present examples of the flare optical LCs exploring different redshifts to the source and compare them with the limiting magnitudes in the  and  bands of ZTF, DECam , and LSST.For DECam search we assumed a 60s exposure 5- detection on a dark night (Bom et al. 2023), while for LSST we assumed the default 30s (2 × 15s) exposures4 (Ivezić et al. 2019).The emission of the emerging cocoon is calculated using a photon diffusion model of a plasma in homologous expansion (see Section 2.2), a model typically employed to interpret LCs of supernovae.Therefore, this EM counterpart evolves qualitatively similar to the emission of supernovae during their most intense phase.However, it can be distinguished from supernovae since the LC of the latter typically lasts from ∼ 40 days to a few months (Wheeler & Harkness 1990;Kasen & Woosley 2009), whereas the emission of the flare discussed here lasts from a few days up to ∼ 40 days, accordingly.
The lag (relative to the GW event) of the flare ranges from weeks to years (see Figures 7 and 8).The flares with the shortest time lags and at the same time magnitude variations |Δ| ≳ 0.5, which we consider the ones that could be better associated to an observed GW event, are produced in AGNs with SMBH ≲ 5 × 10 6 M ⊙ , by remnants with masses ≲ 100 M ⊙ , from mergers occurring at distances  ≳ 4000 g from the central SMBH.

SUMMARY AND DISCUSSION
Active galactic nuclei (AGNs) have been proposed as plausible sites hosting a sizable fraction of the BBH mergers producing the GW events detected by the LIGO-Virgo-Kagra (LVK) experiment (Ford & McKernan 2022).Previous analyses suggest that such GW events could be accompanied by EM counterparts due to the interaction of the merger remnant with the AGN disc gas (Stone et al. 2017;McKernan et al. 2019).In this paper, we work out a new astrophysical scenario leading to thermal flares that result from the interaction of a BBH merger remnant with an AGN thin disc.The proposed scenario is based on the following considerations.We overplot the estimated limiting magnitudes of the ZTF (Bellm et al. 2019), DECam (Bom et al. 2023), and Vera C. Rubin (Ivezić et al. 2019) instruments for the  and  bands (solid and dashed horizontal lines, respectively).The curves in the left and right panels correspond to AGNs with SMBHs of 5 × 10 6 and 5 × 10 7 M ⊙ , respectively, whereas the curves in the upper, middle, and lower panels consider remnants of  BH = 50, 100, and 200 M ⊙ , respectively.In each panel, the upper, middle, and lower curves correspond to sources at redshifts of  = 0.3, 0.6, and 1.3, respectively.We explore LC solutions of remnant-disc interactions occurring at  = 1500 and 4000  g .These solutions correspond to the left and right curve families, at a given redshift.In the panels with one curve family at a given redshift, the LCs correspond to  = 4000 g , since the flare cannot be produced at  = 1500 g within the present analysis (see the condition 4 in Section 2).
• A recoiling and highly spinning BH is formed due to a second or higher generation of BBH merger taking place at a few thousand Schwartzchild radii from the SMBH and outside the disc.
• The kicked BH enters the dense region of the disc, accretes material, and launches relativistic jets that propagate quasi-parallel to the disc plane.
• One of the jet cocoons created within the disc, emerges and expands outside the disc on the observer's side.We then calculate the emission of photons that diffuse and emanate from the surface of the emerging cocoon.
Given the conditions of the cocoon when breaking out the disc, we estimate its emission using a photon diffusion model typically employed to describe LCs of supernovae (Arnett 1996;Chatzopoulos et al. 2012).The emerging cocoon mainly emits at optical and EUV wavelengths and its LC profiles exhibit in general the following features (see Figures 6, 7, and 8): (i) The larger the mass and accretion rate of the SMBH, the larger the time delay (after the GW) and duration of the flare.
(ii) The smaller the remnant's mass, the larger the duration of the flare.
(iii) The larger the kick velocity, the larger the flare duration.
In AGNs with SMBHs of ∼ 10 6−7 M ⊙ , remnants with kick velocities in the range of ∼[300 -1000] Km s −1 can drive flares comparable to or exceeding the emission of the hosting AGN.Such flares exhibit durations and periods of appearance in different time scales, depending on the parameter configuration of the system (see Figures 6 and  7).For instance, in AGNs with SMBHs of ∼ 5 × 10 6 M ⊙ , remnants with masses ≲ 100 M ⊙ interacting with the disc at ∼5000  g from the SMBH can produce flares representing more than 100% of the AGN background emission.These flares are relatively fast, appearing within ∼[10-100] days after the GW and lasting for ∼[1-5] days.Alternatively, remnants with masses ≳ 100 M ⊙ can produce flares exceeding 200% of the background emission when interacting with the disc at distances ≳ 5000 g from SMBHs of ∼ 5×10 7 M ⊙ .These bright flares can appear within [100-1000] days after the GW and last about  days (see Figures 6 and 7).
AGNs typically exhibit stochastic variability in optical bands with |Δ| of a few tenths of magnitudes in timescales from days to years (MacLeod et al. 2012;Graham et al. 2017;Aranzana et al. 2018).Thus, among the possivel flare profiles derived here, we suggest that those with the shortest time lags and magnitude variations |Δ| ≳ 0.5 are the ones that could be better associated with a GW event.Such flares are produced in AGNs with SMBHs of masses ≲ 5 × 10 6 M ⊙ , by remnants with masses ≲ 100 M ⊙ , from mergers occurring at distances ≳ 4000 g from the central SMBH (see Figures 6, 7, and 8).These flares appear within ∼[10-100] days after the GW event, lasting about 5 days, and require remnants with kicks exceeding 600 km s −1 .Flares in AGNs with SMBHs of such size can be detected up to  ∼ 0.25 by ZTF, up to  ∼ 0.9 by DECam, and up to  ∼ 1.3 by LSST (see Figure 8).
The present analysis is appropriate for BBH merger remnants with masses ≳ 50 M ⊙ , consistent with the situation of a second or highergeneration merger in a hierarchical sequence.Motivated by the fact that recoils from previous coalescence can perturb the alignment of the BHs from the disc plane, we assume the BH remnant of the present analysis was born outside the disc.Nevertheless, the present model could also be applied to coalescences occurring within the disc, provided that the pre-merger cavity (Kimura et al. 2021) is much smaller than the disc thickness.
There have been reported so far over seven optical flares candidates associated with BBH GW events in AGNs (Graham et al. 2020(Graham et al. , 2023)).Although none of them are yet confirmed, their theoretical interpretation is timely to favour or disfavour such associations.In a forthcoming work, we will investigate whether the current EM counterparts candidates can be explained in terms of the emission scenario discussed here.
The present emission scenario can also be employed to interpret optical-EUV flares from AGNs with no associated GW event.Such emissions would correspond to flares produced by BHs orbiting around the AGN central engine that eventually threads the disc and emerges on the observer's side.In this case, there would be no observational constraint for the flare starting time, contrary to the case of EM flares with an associated GW signal.
The thermal emission model discussed in this paper relies on an jet cocoon that emerges as a non-relativistic and quasi-spherical outflow of matter outside the AGN disc.Here we propose that such emerging cocoons are more likely driven by jets propagating within and quasiparallel to the disc rather than jets propagating in the perpendicular direction.The latter case would lead to jetted and faster emerging cocoons for which the spherical expansion approximation assumed here would not apply.It remains as subject of further investigation to assess how frequently BBH remnants produce relativistic jets propagating quasi-parallel to the disc.
Clear limitations of the present analytic model are the assumed morphology for the jet cocoon and the disc structure.Given the ejection of relativistic jets within the disc, the jet cocoons can be bent by ram pressure due to the motion of the BH remnant.In addition, the disc vertical density stratification (not considered in the present analysis) is expected to influence the jet morphology and dynamics and hence the properties of the emerging cocoon.Such effects could be accounted with 3D magneto-hydrodynamical simulations of the problem discussed here.Thus, if BBH mergers indeed produce measurable EM counterparts, the predictions of the present analytic model can be employed to benchmark limit cases of more realistic calculations seeking to reconstruct the astrophysical processes leading to multimessenger emission from BBH mergers.set of parameters one obtains through equations (A11-A12)  k = 0.389, 0.338, 0.309, and the minimum escaping radii of 12121, 9143, and 7641  g , respectively.From the lower panel of Figure A1, we note that remnants with  k ∼[100-400] Km s −1 spend among ∼50% and 5% of their orbit time within the disc, whereas remnants with  k ≳ 1500 Km s −1 spend ∼ 5% or less of their time within the disc.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Sketch of the scenario proposed in this work for the EM counterpart to a BBH GW event.The blue rectangular regions in panels (A) and (B) represent the region of the AGN thin disc of half-thickness ℎ d at a distance  from the SMBH, nearby the location of the BBH merger.(A): The BBH coalescence occurs at a vertical distance  0 from the disc mid-plane, leaving a BH remnant with a gravitational kick of velocity  k .The remnant BH launches efficient jets in a period    HL / k after entering the disc.(B): The jets inflate cocoons of shocked disc material (represented by the dark blue regions) that propagate within the disc.(C): The material of one of the jet cocoons emerges and expands outside the disc on the observer side.(D): The emerging cocoon let escape thermal photons producing an observable flare.

Figure 4 .
Figure 4. SEDs of the emerging cocoon emission at different times after the assumed GW event (blue curves).The curves were obtained through the model described in Section 2 considering a merger remnant of 100 M ⊙ with a kick velocity of 800 Km s −1 interacting with a thin disc at a distance of  = 5000 g from the central SMBH.Left and right panels assume an SMBH of 6 × 10 6 and 3 × 10 7 M ⊙ , respectively, whereas the upper and lower panels consider accretion rates of 0.03 and 0.1 (in Eddington units), respectively.In each panel, the dashed black curve plots the emission of the associated thin disc.The grey curves plot the background AGN emission adapted from Ho 2008 (see the text).The source is assumed to be located at a redshift of  = 0.438.For reference, vertical lines indicate examples of frequencies at NIR (J band, 1235 nm, brown), optical (g-band, 464 nm, green), and extreme ultraviolet (UVw2 band, 193 nm, cyan).

Figure 6 .
Figure 6.Fractional excess of the emerging cocoon emission relative to the AGN background (see the text) as a function of the mass of the BH remnant.As indicated, the curves with different line styles correspond to different kick velocities.The curves in brown, green and blue correspond to the fluxes at NIR, optical, and extreme ultraviolet frequencies, respectively, considered in Figure5.The curves in the left panels were calculated with different SMBH masses, as indicated using the fixed values of  = 0.05 and  = 5000 g .The curves in the right panels were calculated with different distances  from the SMBH as indicated with the fixed values of  = 0.05 and  SMBH = 10 7 M ⊙ .

Figure 7 .
Figure 7. Time delay for the appearance of the flare after the GW event (thick curves) and its duration (thin curves) as functions of the mass of the remnant.The curves correspond to different BH kick velocities, as labelled.The curves in the upper panels were obtained using  = 0.05 and  = 5000 g and different masses for the SMBH, as indicated.The curves in the lower panels were calculated considering BBH coalescences at different distances from the SMBH, as indicated, and using the fixed values of  = 0.05 and  SMBH = 10 7 M ⊙ .The source is assumed at redshift  = 0.438.
Time delay Δ pGW (in the observer frame) as a function of the remnant kick velocity for the appearance of the emerging cocoon flare after the GW event.This time interval is a sequence of different periods of the present emission model (see the text), namely the periods in which (i) the BH remnant enters the disc boundary, Δ ent , (ii) the jets form within the disc, Δ j , (iii) the jet cocoon reaches the disc boundary, Δ c , and (iv) the emerging cocoon produces the maximum bolometric luminosity Δ max (see the text).The blue curve plots Δ pGW as the sum of the aforementioned period components.
d , respectively.We see that remnants with kick velocities of ∼ Time delay for the appearance of the flare after the BBH coalescence (measured at Earth) derived as a function of the kick velocity  k of the merger remnant.