Post-merger Jets from Supermassive Black Hole Coalescences as Electromagnetic Counterparts of Gravitational Wave Emission

Supermassive black hole (SMBH) mergers are ubiquitous in the history of the universe (Begelman et al. 1980; Kormendy & Ho 2013; Kroupa et al. 2020) and can produce powerful gravitational wave (GW) bursts when they coalesce (e.g., Thorne & Braginskii 1976; Sesana et al. 2004), making them promising candidates for GW detectors such as the Laser Interferometer Space Antenna (LISA, Amaro-Seoane et al. 2017; Baker et al. 2019) and pulsar timing arrays (PTAs, e.g., Mingarelli et al. 2017; Taylor et al. 2019; Arzoumanian et al. 2020) in single-source and/or stochastic GW background searches. The accretion activity between the binary system and the surrounding disk can produce multiwavelength electromagnetic (EM) emission (e.g., Milosavljević & Phinney 2005; Moesta et al. 2012; Farris et al. 2015; Kelly et al. 2017; Haiman 2017; dAscoli et al. 2018), and the time-variable EM signatures from the circumbinary disks could be detectable (e.g., Schnittman & Krolik 2008; Haiman et al. 2009; Tanaka & Menou 2010). The spinning SMBH expected to form after the SMBHs have coalesced may also lead to relativistic jets, in which particle acceleration will take place. The resulting nonthermal emission from the accelerated electrons may provide a promising post-merger EM counterpart of the GW emission, and will not only provide complementary information on SMBH mergers but also shed light on the physical processes in these systems (e.g., Schnittman 2011; Ravi & Array 2018; Mészáros et al. 2019). Yuan et al. 2020 recently suggested that the SMBH mergers can also be high-energy neutrino emitters, and demonstrated that they are also promising targets for high-energy multi-messenger astrophysics (Murase & Bartos 2019).

We study the EM emission produced in relativistic jets launched after the coalescence of SMBHs. The physical picture is that the disk winds originating from the circumbinary disk and mini-disks around each SMBH form a pre-merger wind bubble, and jets powered by the Blandford–Znajek (BZ; Blandford & Znajek 1977) mechanism are launched after the merger. The jets push ahead inside the pre-merger disk wind material, resulting in the formation of forward and reverse shocks. In the forward shock region, electrons are accelerated to high energies with a power-law distribution as observed in afterglows of gamma-ray bursts (GRBs; e.g., Mészáros 2006). These particles then produce broadband nonthermal EM emission through synchrotron and synchrotron self-Compton (SSC) processes. This Letter is organized as follows. In Section 2, we introduce the physical conditions of the pre-merger wind bubble and model the propagation of jets. The radiation processes and the resulting photon spectra, light curves, and detection horizons are presented in Section 3. In Section 4, we discuss implications of our results. Throughout the Letter, we use the conventional notation Q_x = Q/10^x and physical quantities are written in centimeter-gram-second units, unless otherwise specified.

We discuss here the physical conditions in a pre-merger circumbinary environment and derive relevant quantities that describe the jet propagation. We consider on-axis observers, which is sufficient for the purpose of this work. The emission region is typically expected to be only mildly relativistic on timescales of interest (the corresponding observation time after the jet launch is T ∼ 10⁵–10⁶ s).

Numerical simulations have demonstrated that binary SMBH mergers can produce jet-like emissions driven by the Poynting outflow (e.g., Kelly et al. 2017). We assume that a jet is launched after the coalescence and subsequently propagates in the wind bubble formed by pre-merger disk winds. Figure 1 schematically illustrates the configuration of the system. The gaseous environment in the host galaxy center, in which the wind bubble expands. We focus on shock interactions between the jet and wind bubble. Initially, the circumbinary disk can react promptly to the evolution of the binary system. The ratio between the disk radius R_d and the semimajor axis of the binary system a remains unchanged (R_d /a ∼ 2), until the inspiral timescale t_GW of the binary system (e.g., Shapiro & Teukolsky 1983) equals the viscosity timescale t_vis (e.g., Pringle 1981), which is known as the disk decoupling. After the disk becomes decoupled, the merger of SMBHs in binary system occurs within the time interval ${t}_{m}\sim 3\,\times {10}^{-2}\,\mathrm{yr}\,{M}_{\mathrm{BH},6}{\alpha }_{-1}^{-8/5}{h}_{-1}^{-16/5}$, where M_BH = 10⁶ M_BH,6 M_⊙ is the mass of the binary system, the dimensionless parameter h is defined as h = H/R_d , α ∼ 0.1 is the viscosity parameter, and H is the disk scale height. The disk gas starts to fill the cavity between the disk and the SMBHs in the viscosity timescale,

$$\begin{eqnarray}&&{t}_{\mathrm{vis}}\sim 0.1\,\mathrm{yr}\,{M}_{\mathrm{BH},6}{\alpha }_{-1}^{-8/5}{h}_{-1}^{-16/5}\end{eqnarray} \tag{ 1 }$$

after the coalescence (Farris et al. 2015). This leads to a time delay (t_delay ∼ t_vis) of days to months between the GW burst and the launch of post-merger jets, if h ∼ 0.1–0.3 is assumed. However, for a thick and highly magnetized disk with h ∼ α ∼ 1, t_delay could be much shorter.

Zoom In Zoom Out Reset image size

On the other hand, within the duration of these two short-term processes, e.g., t_m and t_vis, the disk wind radius may reach v_d (t_vis − t_m ) ∼ 10¹⁴ − 10¹⁶cm above the disk, where v_d is the disk wind velocity that is of the order of the escape velocity, ${v}_{\mathrm{esc}}({R}_{d,\mathrm{dec}})\approx \sqrt{2{{GM}}_{\mathrm{BH}}/{R}_{d,\mathrm{dec}}}$ for the circumbinary disk, and ${R}_{d,\mathrm{dec}}\approx 1.2\times {10}^{13}\,\mathrm{cm}\,{M}_{\mathrm{BH},6}{\alpha }_{-1}^{-2/5}{h}_{-1}^{-4/5}$ is the radius of the circumbinary disk at the decoupling. In reality, not only the circumbinary disk but also mini-disks around two SMBHs contribute, which would make the wind bubble more complicated. For simplicity, we assume the density profile of the winds at the decoupling to obtain the density distribution of the wind bubble at larger distances,

$$\begin{eqnarray}&&{\varrho }_{w}(r)=\displaystyle \frac{{\eta }_{w}(1+\chi ){\dot{M}}_{\mathrm{BH}}}{4\pi {r}^{2}{v}_{d}}\equiv {{Dr}}^{-2},\end{eqnarray} \tag{ 2 }$$

where ${\dot{M}}_{\mathrm{BH}}$ is the mass accretion rate onto the binary system, χ ∼ 1 is introduced to take into account the contribution of mini-disks, and η_w represents the fraction of accreted mass converted to the disk wind. According to the simulations, for SANE (Standard And Normal Evolution) models, the parameter η_w may vary from 10⁻⁴ to 10⁻¹ (Jiang et al. 2019b, 2019a; Ohsuga et al. 2009) when the mass accretion rate changes from sub-Eddington to super-Eddington. In MAD (Magnetically Arrested Disk) models, η_w can reach 10⁻² to 10⁻¹ (Akiyama et al. 2019). With v_d ∼ v_esc(R_d,dec), we have $D\simeq 5.9\times {10}^{11}\,{\rm{g}}\,{\mathrm{cm}}^{-1}\,{\tilde{\eta }}_{w,-1.5}(\dot{m}/0.5){M}_{\mathrm{BH},6}{\beta }_{d,-1}^{-1}$, where ${\tilde{\eta }}_{w}\,\equiv (1+\chi ){\eta }_{w}$, β_d,−1 ≡ v_d /(0.1c), the parameter $\dot{m}$ is defined as the ratio of ${\dot{M}}_{\mathrm{BH}}$, and the Eddington value ${\dot{M}}_{\mathrm{Edd}}\,\equiv 10{L}_{\mathrm{Edd}}/{c}^{2}$ (assuming a radiation efficiency of 0.1).

After the coalescence, a powerful jet driven by the spin energy of the newly formed SMBH can appear, subsequently propagating in the pre-merger wind bubble. Considering a sub-Eddington accretion rate with the MAD configuration, we estimate the jet kinetic luminosity to be

$$\begin{eqnarray}\begin{array}{rcl}{L}_{k,j} & = & {\eta }_{j}{\dot{M}}_{\mathrm{BH}}{c}^{2}\\ & \simeq & 6.3\times {10}^{44}\,\mathrm{erg}\ {{\rm{s}}}^{-1}{\eta }_{j}(\dot{m}/0.5){M}_{\mathrm{BH},6},\end{array}\end{eqnarray} \tag{ 3 }$$

where η_j ∼ 0.3–1 is the ratio of the accretion energy converted to the jet energy (Tchekhovskoy et al. 2011). Following the standard jet propagation theory (e.g., Bromberg et al. 2011; Mizuta & Ioka 2013), we write down the dimensionless parameter that represents the ratio of the energy density of the jet and the rest-mass energy density of the surrounding medium

$$\begin{eqnarray}&&\tilde{L}\approx \displaystyle \frac{{L}_{k,\mathrm{iso}}}{4\pi {r}^{2}{\varrho }_{w}{c}^{3}}\simeq 63\,{\tilde{\eta }}_{w,-1.5}^{-1}{\eta }_{j}{\theta }_{j,-0.5}^{-2}{\beta }_{d,-1},\end{eqnarray} \tag{ 4 }$$

where θ_j is the jet opening angle, ${L}_{k,\mathrm{iso}}\approx 2{L}_{k,j}/{\theta }_{j}^{2}$ is the isotropic-equivalent luminosity. Since the quantity $\tilde{L}$ lies in the regime ${\theta }_{j}^{-4/3}\ll \tilde{L}$, we expect that the jet is "uncollimated" for our fiducial parameters. ⁶ This situation is similar to that in choked jet propagation in the circumstellar material (Senno et al. 2016; Nakar 2015), and β_h is evaluated from $\tilde{L}$. In the relativistic limit, the jet head Lorentz factor is ${{\rm{\Gamma }}}_{h}\approx {\tilde{L}}^{1/4}/\sqrt{2}$ (Senno et al. 2016), and we have Γ_h ∼ 2 in our fiducial case with ${\tilde{\eta }}_{w}={10}^{-1.5}$. Note that the jet head radius is ${R}_{h}=c{\beta }_{h}\hat{T}\approx c\hat{T}$, and $\hat{T}$ is introduced to represent time measured in the central engine frame, which can be converted to the observation time T via $T=(1+z)(1-{\beta }_{h})\hat{T}$ (that is $T\approx (1+z)\hat{T}/[2{{\rm{\Gamma }}}_{h}^{2}]$ in the relativistic limit) for on-axis observers.

Furthermore, to ensure particle acceleration, we impose radiation constraints requiring that the shock is collisionless, without being mediated by radiation (Murase & Ioka 2013; Senno et al. 2016). Here, ignoring effects of pair production, we use the conservative condition, τ_T ≈ ϱ_w σ_T R_h /m_p < 1, where σ_T is the Thomson cross section. Numerically, this condition is satisfied at $\hat{T}\gtrsim 10\,{\rm{s}}$, which is much shorter than the duration of EM emission.

With the jet dynamics presented in the previous section, we calculate the EM spectra resulting from synchrotron and SSC emission. As in the standard theory of GRB afterglows (e.g., Mészáros 2006), we assume that electrons are accelerated at the external forward shock with a power-law spectral index s. The energy fractions of the downstream energy density converted to nonthermal electron and magnetic field energy are defined as _e and _B , respectively. The upstream number density is given by ${n}_{h,u}\approx {\varrho }_{w}({R}_{h})/{m}_{p}\propto {R}_{h}^{-2}$, and $B\,\approx {[{\epsilon }_{B}32\pi {{\rm{\Gamma }}}_{h}({{\rm{\Gamma }}}_{h}-1){n}_{h,u}{m}_{p}{c}^{2}]}^{1/2}$ is the downstream magnetic field strength.

In the relativistic limit (Γ_h ≫ 1), the characteristic injection frequency ν_m and the cooling frequency ν_c in the observer frame are written, respectively, as,

$$\begin{eqnarray}\begin{array}{rcl}{\nu }_{m} & \approx & \displaystyle \frac{3{{\rm{\Gamma }}}_{h}{\gamma }_{m}^{2}{eB}}{4\pi (1+z){m}_{e}c}\\ & \simeq & 3.4\times {10}^{3}\,\mathrm{GHz}\ {\epsilon }_{e,-1}^{2}{\zeta }_{e,-0.4}^{2}{\epsilon }_{B,-2}^{1/2}{\eta }_{j}^{1/2}{T}_{4}^{-1}\\ & & \times {\left(\dot{m}/0.5\right)}^{1/2}{M}_{\mathrm{BH},6}^{1/2}{\theta }_{j,-0.5}^{-1}\end{array}\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{\nu }_{c} & \approx & \displaystyle \frac{3{{\rm{\Gamma }}}_{h}{\gamma }_{c}^{2}{eB}}{4\pi (1+z){m}_{e}c}\\ & \simeq & 4.6\times {10}^{2}\,\mathrm{GHz}\ {\left(1+z\right)}^{-2}{\left(1+Y\right)}^{-2}{\epsilon }_{B,-2}^{-3/2}\\ & & \times {\tilde{\eta }}_{w,-1.5}^{-2}{\theta }_{j,-0.5}^{-1}{\left(\dot{m}/0.5\right)}^{-3/2}{T}_{4}{\beta }_{d,-1}^{2}{M}_{\mathrm{BH},6}^{-3/2},\end{array}\end{eqnarray} \tag{ 6 }$$

where γ_m = _e ζ_e (Γ_h − 1)m_p /m_e is the electron minimum Lorentz factor, and ${\gamma }_{c}=6\pi {m}_{e}c/[(1+Y){T}^{{\prime} }{\sigma }_{T}{B}^{2}]$ is the cooling Lorentz factor. Here, ${\zeta }_{e}={g}_{s}/{f}_{e}=1/[{f}_{e}\mathrm{ln}({\gamma }_{M}/{\gamma }_{m})]\,\sim 0.3\mbox{--}0.4$ is constrained by the particle-in-cell simulations (Park et al. 2015) (where f_e is the fraction of accelerated electrons and the maximum Lorentz factor of electrons is ${\gamma }_{M}={(6\pi e)}^{1/2}/{[{\sigma }_{T}B(1+Y)]}^{1/2}$), Y is the Compton parameter, and ${T}^{{\prime} }=\hat{T}/{{\rm{\Gamma }}}_{h}\approx 2{{\rm{\Gamma }}}_{h}T/(1+z)$ is the comoving time. For example, at T = 10⁴ s, we have Y ≃ 2.4, corresponding to the fast cooling regime. It changes to the slow cooling regime on a timescale from days to weeks. We obtain the peak synchrotron flux (e.g., Wijers & Galama 1999)

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\nu ,\mathrm{syn}}^{\max } & \approx & \displaystyle \frac{(1+z)(0.6{f}_{e}{n}_{h,u}{R}_{h}^{3}){{\rm{\Gamma }}}_{h}{e}^{3}B}{\sqrt{3}{m}_{e}{c}^{2}{d}_{L}^{2}}\\ & \simeq & 0.24\,{\rm{m}}\ \mathrm{Jy}\,(1+z){g}_{s,-1.2}{\zeta }_{e,-0.4}^{-1}{\left(\dot{m}/0.5\right)}^{3/2}\\ & & \times {\eta }_{j}^{1/2}{\epsilon }_{B,-2}^{1/2}{\tilde{\eta }}_{w,-1.5}{\beta }_{d,-1}^{-1}{\theta }_{j,-0.5}^{-1}{M}_{\mathrm{BH},6}^{3/2}{d}_{L,28}^{-2}.\end{array}\end{eqnarray} \tag{ 7 }$$

The low-frequency synchrotron emission is subject to synchrotron self-absorption (SSA). The SSA optical depth is written as ${\tau }_{\mathrm{ssa}}{(\nu )={\xi }_{s}{{en}}_{h,u}{R}_{h}(\nu /{\nu }_{n})}^{-p}/[B{\gamma }_{n}^{5}]$, where ν is the observed frequency, ξ_s ∼ 5–10 depends on the electron spectral index, ${\gamma }_{n}=\min \left[{\gamma }_{m},{\gamma }_{c}\right]$, ${\nu }_{n}={\gamma }_{n}^{2}{eB}/[(1+z){m}_{e}c]$, p = 5/3 for ν < ν_n and p = (4 + s)/2 or p = 3 for ν > ν_n depending on the slow or fast cooling regime (e.g., Panaitescu & Kumar 2000; Murase et al. 2014). The critical timescales set by τ_ssa = 1 for ν < ν_n and for ν > ν_n are ${T}_{\mathrm{ssa}}\simeq 5.4\times {10}^{5}\,{\rm{s}}\,{\xi }_{s,1}^{3/10}{(1+z)}^{1/2}{(1+Y)}^{1/2}\,{\epsilon }_{B,-2}^{3/5}{\left(\tfrac{\nu }{1\mathrm{GHz}}\right)}^{-1/2}{\tilde{\eta }}_{w,-1.5}^{11/10}{\beta }_{d,-1}^{-11/10}{(\dot{m}/0.5)}^{9/10}{M}_{\mathrm{BH},6}^{9/10}$ and ${T}_{\mathrm{ssa}}\simeq 3.5\times {10}^{5}\,{\rm{s}}\,{\xi }_{s,1}^{1/2}{(1+z)}^{-1/2}{(1+Y)}^{-1/2}{M}_{\mathrm{BH},6}^{1/2}{\left(\tfrac{\nu }{100\mathrm{GHz}}\right)}^{-3/2}{(\dot{m}/0.5)}^{1/2}{M}_{\mathrm{BH},6}^{1/2}{\tilde{\eta }}_{w,-1.5}^{1/2}{\beta }_{d,-1}^{-1/2}$, respectively. Thus, we expect that EM emission at 5 GHz and 100 GHz reaches a peak about a few days after the jet launch (T_ssa ≃ 7.1 × 10⁵ s and T_ssa ≃ 3.1 × 10⁵ s, respectively, in our fiducial case with ξ_s = 8.7).

We numerically calculate the electron distribution and the resulting synchrotron and SSC spectra of the forward shock, following the method used in Murase et al. (2011) and Zhang et al. (2020). We solve the continuity equation that describes the evolution of the electron spectra and calculate the synchrotron/SSC components, in which the trans-relativistic regime can be consistently treated as in Zhang et al. (2020). Combining the obtained radio, millimeter, optical and X-ray light curves with the sensitivities of corresponding detectors, we discuss the possibility of follow-up observations of the EM counterpart.

The left panel of Figure 2 shows the snapshots of synchrotron and SSC spectra at T = 10⁴ − 10⁶ s for an on-axis source located at z = 1. We assume s = 2.0, _e = 0.1, and _B = 0.01. The solid and dashed lines correspond to the synchrotron and SSC components. Very high-energy gamma-ray emission at ≳1 TeV energies is suppressed due to the Klein-Nishina effect (e.g., Murase et al. 2011; Zhang et al. 2020), and the γ γ annihilation with the extragalactic background light (EBL). For the EBL correction, γ γ optical depth data from Model C in Finke et al. (2010) is used. To show how the EM signal evolves with time, we illustrate the gamma-ray (1 GeV), X-ray (1 keV), UV (1 eV), and radio (5 GHz and 100 GHz) light curves in the right panel. In particular, before the characteristic time T_ssa (shown as the vertical yellow and blue lines). The radio emission is suppressed by the SSA process, which is implemented by multiplying $(1-{e}^{-{\tau }_{\mathrm{ssa}}})/{\tau }_{\mathrm{ssa}}$.

Zoom In Zoom Out Reset image size

It is useful to discuss the detection horizon ${d}_{\mathrm{lim}}$ for some detectors such as the Square Kilometre Array (SKA), Very Large Array (VLA), Expanded VLA (EVLA), Atacama Large Millimeter Array (ALMA), Hubble Space Telescope (HST), James Webb Space Telescope (JWST), Large Synoptic Survey Telescope (LSST), and the high-resolution camera on the Chandra X-ray Observatory (Chandra) ⁷ as functions of the observation time T. Given the observed flux F_ν (ν_γ , T, z) at the observer time T from an on-axis source located at redshift z, the horizon can be calculated iteratively via

$$\begin{eqnarray}&&{d}_{\mathrm{lim}}({\nu }_{\gamma },T)={d}_{L}{\left(\displaystyle \frac{\tfrac{1}{{\rm{\Delta }}{T}_{\exp }}{\int }_{T}^{T+{\rm{\Delta }}{T}_{\exp }}{F}_{\nu }({\nu }_{\gamma },t,z){dt}}{{F}_{\mathrm{lim}}({\nu }_{\gamma },{\rm{\Delta }}{T}_{\exp })}\right)}^{1/2},\end{eqnarray} \tag{ 8 }$$

where ${F}_{\mathrm{lim}}({\nu }_{\gamma },{\rm{\Delta }}{T}_{\exp })$ is the detector sensitivity normalized to the exposure time ${\rm{\Delta }}{T}_{\exp }$. For example, specifying the detection frequency ν = 100 GHz, the sensitivity of ALMA is approximately 34 μJy for a 1 hr integration, e.g., ${\rm{\Delta }}{T}_{\exp }=1$ hr. Figure 3 indicates the detection horizons for SKA (5 GHz, ${\rm{\Delta }}{T}_{\exp }=10$ hr), SKA (1 GHz ⁸ , ${\rm{\Delta }}{T}_{\exp }=10$ hr), VLA (5 GHz, ${\rm{\Delta }}{T}_{\exp }=1$ hr), ALMA (100 GHz, ${\rm{\Delta }}{T}_{\exp }=1$ hr), JWST (1 eV, ${\rm{\Delta }}{T}_{\exp }=10$ ks), HST (1 eV, ${\rm{\Delta }}{T}_{\exp }=10$ ks), LSST (r-band, point source exposure time ${\rm{\Delta }}{T}_{\exp }=30\,{\rm{s}}$ in the 3 day revisit time), and Chandra (1 keV, ${\rm{\Delta }}{T}_{\exp }=100$ ks). The vertical black and blue dotted lines, respectively, illustrate the times T_ssa at which photons at 100 GHz and 5 GHz bands start to survive from the synchrotron self-absorption.

Zoom In Zoom Out Reset image size

From Figure 3, we expect that ALMA, SKA, and EVLA can detect SMBH mergers in the radio bands, respectively, out to redshifts of z ∼ 4–6. Remarkably, the optical and X-ray signals from the mergers in the range of 1 ≲ z ≲ 2 can also be identified through targeted searches by Chandra, HST, and JWST in a long duration after the merger. In addition, we can estimate the observation time for each detector if the luminosity distance of the merger is specified. For example, a source located at z = 3 would remain detectable by ALMA for roughly 20–30 days (see the black dotted horizontal line). One caveat is that this calculation is carried out in the ideal case where the detectors can point to the position of the source and start the observation immediately after the EM signal reaches the Earth. We discuss the sky coverage and a detection strategy in the following Section 4.

We investigated broadband nonthermal EM emission from electrons accelerated at the external forward shock expected in post-merger jets from the coalescence of SMBHs. In our model, the jets can be launched at t_delay ∼ t_vis ∼ (0.003–0.1)M_BH,6 yr after the coalescence. The time lag is primarily determined by the scale height of the circumbinary disk and the viscosity parameter. We found that, for a moderate accretion rate ($\dot{m}\sim 0.5$), the multiwavelength emission from such a system may persist at detectable levels for months after the jet launch, depending on the facilities and the luminosity distance. Moreover, according to our model, the sources with moderate $\dot{m}=0.5$ can be detected up to z ∼ 5–6, covering the range that LISA-like GW detectors have the best detection chance, e.g., z ∼ 1–2, in which (1–10)f_b mergers per year are expected (Menou et al. 2001; Enoki et al. 2004; Arun et al. 2009; Amaro-Seoane et al. 2012; Dal Canton et al. 2019; Kara et al. 2019). Here ${f}_{b}\sim 1/(2{{\rm{\Gamma }}}_{h}^{2})$ is the beaming factor in our model. Because the jet head Lorentz factor is as low as Γ_h ≲ 2, the EM emission from the forward shock region is not highly beamed and we expect f_b ∼ 0.1–1. This makes the binary SMBH mergers interesting targets for future multi-messenger studies. If super-Eddington accretion (e.g., $\dot{m}\sim 10$) occurs, as was optimistically assumed in Yuan et al. (2020), even LSST and Chandra could detect EM signals from the sources in the redshift range 4 ≲ z ≲ 6. We showed the case of s = 2.0 for the demonstration. If a larger spectral index, e.g., s ∼ 2.2–2.4, is used, as expected from observations of GRB afterglows, the radio detection would be more promising whereas a higher accretion rate would be required for successful optical and X-ray observations.

The density of the pre-merger bubble, which was assumed to be a wind profile, is subject to large uncertainties. The extrapolation in the density distribution would be applicable up to an outer wind radius of ∼10¹⁴–10¹⁶ cm. The density predicted by Equation (2) would drop below that of the central molecular zone (indicated as the circumnuclear environment in Figure 1), which may lead to an increased radio emission. In addition, a cocoon formed along with the jet, depending on uncertain details of the medium, could produce thermal photons, which may not only lead to detectable signals but also serve as seed photons for inverse-Compton emission. We focused on the more secure EM emission from the forward shock region as the jet propagates in the wind. In this sense, our prediction for the fluxes are conservative.

EM emission from the external reverse shock and internal shocks can also be expected (e.g., Mészáros & Rees 1999; Kobayashi & Zhang 2003, for the reverse shock emission in GRBs). Qualitatively, the ratio between the peak fluxes of the reverse and forward shock emission depends on the value of Γ_j , and the reverse shock contribution might be important for Γ_j ≫ Γ_h .

Previous studies based on general relativistic three-dimensional magnetohydrodynamics simulations have shown that the circumbinary disk and the corona can emit light in UV/EUV bands (e.g., dAscoli et al. 2018), while X-ray and infrared emission from the post-merger circumbinary disk are expected to last for years (Milosavljević & Phinney 2005; Schnittman & Krolik 2008). In the pre-merger phase, the orbits of dual SMBH cores may be identified by radio facilities such as VLBI (e.g., Rodriguez et al. 2006). Blind searches could identify radio or UV/EUV sources from the binary SMBH systems, which would provide complementary constraints on the source location, the accretion rate, and the ambient gaseous environment.

Our model can provide a guidance, including the onset times and the detection windows, in developing detection strategies for future EM follow-up observations, once GW signals are detected. Considering the large uncertainties in the localization with GW detectors, an initial follow-up using large field-of-view (FOV) telescopes, like SKA and LSST, would be necessary to more precisely localize the position of the source. After that, we can use the putative positional information from the initial follow-up imaging to guide the observation of narrower FOV telescopes. In particular, for high-redshift mergers in the range z ∼ 2–5, EM follow-up observations rely more on radio detectors, and the detection is possible a few weeks after the merger. SKA needs the source localization before follow-up observations by VLA and ALMA. On the other hand, if the merger is close enough (e.g., z ∼ 1), LISA observations starting from a few weeks before the merger can localize the merger with a median precision of ∼1 deg² (Mangiagli et al. 2020). In this case, LISA and LSST can jointly guide other X-ray and optical facilities in the very early stage. Amid these two regimes, e.g., z ∼ 1–2, detections in the optical and X-ray bands using HST, JWST, and Chandra would be promising once the source is localized by SKA.

C.C.Y. and P.M. acknowledge support from the Eberly Foundation. The work of K.M. is supported by NSF grant No. AST-1908689, and KAKENHI No. 20H01901 and No. 20H05852. B.T.Z. acknowledges the IGC fellowship. S.S.K. acknowledges the JSPS Research Fellowship, JSPS KAKENHI grant No. 19J00198.