© The Authors 2023

1. Introduction

The merger of two neutron stars detected as both a gravitational wave (GW) and electromagnetic (EM) source (Abbott et al. 2017) has recently opened up the field of multi-messenger studies of astrophysical phenomena. Another promising candidate for such multi-messenger studies is the merger of two massive black holes (BHs). Merging BHs with masses ∼10⁴ − 10⁷ M_⊙ can be detected in GWs by the future space-based interferometer LISA (Laser Interferometer Space Antenna, Amaro-Seoane et al. 2023) as well as by the proposed missions TianQin (Luo et al. 2016) and Taiji (Ruan et al. 2020a). The horizon of these detectors is large, with the ability to detect merging BHs out to z ∼ 10 for mass ratios not too far from unity. Mergers of such massive BHs are also expected to be detectable electromagnetically, as massive BHs are generically surrounded by gas in galactic centres and they are therefore associated with luminous sources such as active galactic nuclei (AGN) when they accrete such gas. If these merging BHs can also be detected electromagnetically, then we can use them to study accretion physics in dynamical spacetimes, to obtain independent measures of BH masses (Amaro-Seoane et al. 2023), and to constrain fundamental physics (Arun 2022) and cosmology (Auclair et al. 2022).

While massive BHs have been detected electromagnetically for many years, in the form of AGN, it remains very unclear whether they give off a sufficiently distinguishable signal at the moment of merger to allow for a multi-messenger study. Over the years, many models and simulations have been developed for the actual physics of the production of an EM counterpart at the merger (e.g., Armitage & Natarajan 2002; Schnittman & Krolik 2008; Rossi et al. 2010; Sesana et al. 2012; Roedig et al. 2014; Gutiérrez et al. 2022; Kelly et al. 2021).

Numerical studies have shown that a change in the luminosity occurs around the time of the merger for BH binaries evolving in circumbinary discs. Before merging, during the late inspiral of gas-rich BH binaries, the binary torques excavate a low-density cavity in the circumbinary disc. Consequently, the circumbinary disc acquires an inner rim at a radius of the order of the BH binary semi-major axis. Material pile sup at this inner rim, creating a high-density region or a non-axisymmetric ‘lump’ (e.g., Kocsis et al. 2012; Noble et al. 2012). Despite the potential barrier that maintains the cavity, gas streams can flow through it and feed accretion minidiscs around the individual BHs (e.g., Noble et al. 2012; Shi & Krolik 2015; Tang et al. 2018). At some point during the binary evolution, the rate at which the orbit shrinks due to the emission of gravitational waves is faster than the viscous timescale – the disc cannot evolve fast enough and decouples from the binary (e.g., Milosavljević & Phinney 2005). Although the gas accretion streams can continue to feed the BHs virtually until the merger, the accretion rate can decrease after the binary completely decouples from the disc (e.g., Gold et al. 2014; Farris et al. 2015). The minidiscs, which dominate the hard X-ray emission (e.g., d’Ascoli et al. 2018), could also gradually disappear near the merger, causing a drop in the X-ray luminosity (Tang et al. 2018). In contrast to this, Armitage & Natarajan (2002) and Cerioli et al. (2016) suggest that squeezing and rapid accretion in the minidisc of the primary BH could lead to large enhancements in the accretion rate and luminosity.

After the BH merger, the disc is expected to maintain initially a central cavity at ∼10 − 20 gravitational radii, causing a small diminution (less than a factor of ∼3 for a non-spinning BH) in the radiative efficiency compared to that of a disc around a single BH (Bogdanović et al. 2022). The cavity is then refilled on a timescale t_vis ∼ 0.1(M_•/10⁶ M_⊙)(α/0.1)^−8/5(h/0.1)^−16/5 yr, where α is the disc viscosity parameter and h is the disc aspect ratio (e.g., Milosavljević & Phinney 2005; Farris et al. 2015; Yuan et al. 2021). The disc refilling gradually increases the BH accretion rate and luminosity. Moreover, when the pileup of material at the inner rim is accreted, the larger availability of gas can potentially drive a post-merger luminosity burst.

Changes in jet properties have also been suggested in conjunction with BH mergers (Merritt & Ekers 2002). The merger can modify BH properties, such as M_•, f_Edd, or a, or the properties of gas dynamics and magnetic fields around the remnant BH, which can lead to observable radio signatures on shorter scales, of the order of hours or days. This is supported by results from general relativistic magnetohydrodynamical simulations of BH mergers (Palenzuela et al. 2010; Moesta et al. 2012; Gold et al. 2014; Kelly et al. 2017, 2021; Cattorini et al. 2021, 2022). For example, simulations predict the production of a flare of Poynting luminosity at the merger, lasting for ∼0.1 day (M_•/10⁶ M_⊙). This Poynting luminosity can be comparable to that of the pre-merger jet, although it seems to depend on various parameters (mass ratio, spin magnitude and alignment, gas density, magnetic field, etc.) which have not been extensively explored in simulations. Yuan et al. (2021) modelled the jet spectrum under the assumption that a newly formed jet after the merger impacts and shocks the nearby gas, thus producing a source similar to a gamma-ray burst.

Fewer studies have been applied to BH merger populations to estimate the number and properties of BH mergers with an EM counterpart (Dotti et al. 2006; Tamanini et al. 2016; Kelley et al. 2019; Krolik et al. 2019; Mangiagli et al. 2022; Lops et al. 2023; Chakraborty et al. 2023). In the following, we study the properties of the BH merger population in the OBELISK simulation (Trebitsch et al. 2021) and assess the multi-messenger observability of their corresponding GW events and EM counterparts, as well as the biases of the observable population. OBELISK is a cosmological radiation-hydrodynamical simulation evolving a protocluster down to redshift ∼3.5. This simulation is ideal for our purposes since it has a high resolution (∼35 pc) and incorporates detailed models for a wide range of BH physical processes, such as accretion, feedback, spin evolution, and dynamical friction, which are key in order to produce a realistic BH merger population. We remind the reader that OBELISK models the evolution of an overdense region, and thus it cannot be used to predict merger rates in an unbiased way. This work follows on from Dong-Páez et al. (2023) in which we present and analyse the population of BH mergers in comparison to the total population of BH in OBELISK.

In Sect. 2, we summarise the properties of the OBELISK simulation and the identification and selection criteria of galaxies and BH mergers. We also describe our calculation of sub-grid merger delays, BH luminosities in several EM bands, and our simulations of the GW parameter estimation by LISA. We present our results in the subsequent sections – in Sect. 3.1, we study the GW observability and parameter estimation by LISA of the BH merger sample, and in Sects. 3.2 and 3.3 we study their observability in several EM bands (X-rays, UV, and radio), the bias in the properties of the observable mergers with respect to the unobserved sample and the synergies with the GW detections. In Sect. 4, we discuss our methods and results in the context of the previous work. Finally, in Sect. 5, we conclude and summarise our main results.

2. Method

2.1. The OBELISK simulation

OBELISK (Trebitsch et al. 2021) re-simulates at high-resolution (Δx ≃ 35pc) the most massive halo in HORIZON-AGN at z = 1.97 in the HORIZON-AGN (Dubois et al. 2014a) volume until redshift z ≃ 3.5. In Fig. 1, we show the projected gas density in a region of the simulation at z ∼ 4. Below we present a brief summary of the properties of the simulation. For a more detailed description, we refer the reader to Trebitsch et al. (2021) and Dong-Páez et al. (2023).

Fig. 1. Gas density projection of a region of the OBELISK simulation at redshift 4.02. BHs with M• > 106 M⊙ are denoted by black star symbols.

The simulation assumes a Λ CDM cosmology with WMAP-7 parameters (Komatsu et al. 2011) – Hubble constant H₀ = 70.4 km s⁻¹ Mpc⁻¹, dark energy density parameter #x03A9;_Λ = 0.728, total matter density parameter #x03A9;_m = 0.272, baryon density parameter #x03A9;_b = 0.0455, amplitude of the power spectrum σ₈ = 0.81, and spectral index n_s = 0.967. The zoomed-in region, with a volume of ∼10⁴ h⁻³ cMpc³, was simulated with a DM mass resolution of 1.2 × 10⁶ M_⊙, while the remaining volume of the original 100 h⁻¹ cMpc HORIZON-AGN box maintained a lower resolution.

OBELISK was run with RAMSES-RT (Rosdahl et al. 2013; Rosdahl & Teyssier 2015), a radiative transfer hydrodynamical code which builds on the adaptive mesh refinement (AMR) RAMSES code (Teyssier 2002). Cells are refined up to a smallest size of 35 pc if its mass exceeds 8 times the mass resolution. The simulation assumes an ideal monoatomic gas with adiabatic index γ = 5/3 and includes gas cooling and heating down to very low temperatures (50 K) with non-equilibrium thermo-chemistry for hydrogen and helium, and contribution to cooling from metals (at equilibrium with a standard ultraviolet background) released by SNe.

Stellar particles have a mass of ∼10⁴ M_⊙, and assume a Kroupa (2001) initial mass function between 0.1 and 100 M_⊙. Stars form in gas cells with density higher than 5 H cm⁻³ and Mach number ℳ ≥ 2. The star formation efficiency depends on the local gas density, sound speed, and turbulent velocity. SN feedback takes place 5 Myr after the birth of a stellar particle, with a mass fraction of 0.2, and releasing 10⁵¹ erg per SN. OBELISK also includes modelling of dust as a passive variable.

BHs form when in a given cell both gas and stars exceed a density threshold of 100 H cm⁻³ and the gas is Jeans unstable. The initial mass is M_•, 0 = 3 × 10⁴ M_⊙ and an exclusion radius of 50 comoving kpc is enforced to avoid formation of multiple BHs. Gas accretion is modelled using Bondi-Hoyle-Lyttleton (BHL) formalism,

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{BHL}}=\frac{4\pi G^2{M}_{\bullet }^2\overline{\rho }}{{({\overline{c_{\mathrm{s}}}}^2+{\overline{v}}_{\mathrm{rel}}^2)}^{3/2}},\end{array}$$(1)

where $\overline{\rho }$, $\overline{c_{\mathrm{s}}}$ and ${\overline{v}}_{\mathrm{rel}}$ are the local average gas density, gas sound speed, and BH relative velocity with respect to the gas. The accretion rate is capped at the Eddington rate

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{Edd}}=\frac{4\pi G{M}_{\bullet }{m}_{\mathrm{p}}}{{\epsilon }_{\mathrm{r}}{\sigma }_{\mathrm{T}}c},\end{array}$$(2)

where m_p is the proton mass, ε_r is the radiative efficiency, σ_T is the Thompson cross-section, and c is the speed of light. A fraction ε_r of the accretion power Ṁc² is radiated, while the remaining 1 − ε_r is accreted onto the BH, contributing to its mass growth. BHs with an Eddington ratio as f_Edd = Ṁ/Ṁ_Edd < 0.01 (here Ṁ = min{Ṁ_BHL,Ṁ_Edd}), are assumed to be radiatively inefficient, and radiative efficiency is reduced by a factor f_Edd/0.01.

AGN feedback is modelled with a dual-mode approach. At f_Edd < 0.01, the AGN releases a fraction ε_MCAF of the rest-mass accreted energy as kinetic energy in jets. Jets assume Magnetically Chocked Accretion Flows, and ε_MCAF is a polynomial fit to the simulations of McKinney et al. (2012) as a function of BH spin. At higher f_Edd ≥ 0.01, the 15% of the accretion luminosity is released isotropically as thermal energy.

Two BHs are merged when their separation becomes smaller than 4Δx and they are gravitationally bound. The simulation models dynamical friction explicitly, including both gas and collisionless particles (stars and DM; following the implementation presented in Dubois et al. 2013; Pfister et al. 2019).

BH spins are self-consistently evolved on the fly via gas accretion and BH-BH mergers following Dubois et al. (2014b). The model for BHs with f_Edd > 0.01 includes evolution of both spin magnitude (Bardeen 1970) and direction (King et al. 2005), while for f_Edd < 0.01) rotational energy is assumed to power the radio jets and therefore the magnitude of BH spins can only decrease. We adopt the polynomial fits in McKinney et al. (2012), with the same procedure for the update of the spin direction as for the f_Edd ≥ 0.01 case. Spin also evolves following the coalescence of two BHs using an analytical fit from Rezzolla et al. (2008). The value of the spin is used to determine the efficiency of the energy injection into jets in the radio mode, and the BH radiative efficiency.

2.2. Galaxy catalogues and BH-galaxy matching

The galaxy and BH merger catalogues used here are identical to those presented in Dong-Páez et al. (2023). We summarise the relevant procedure here and in the next section, but refer the reader to that paper for further details.

Galaxies and their DM haloes were identified together, using a version of the ADAPTAHOP algorithm (Aubert et al. 2004; Tweed et al. 2009) designed to work on both stars and DM particles. Substructures were identified using the most massive sub-maximum method, with a minimum number of particles (stars + DM) of 100 (see details in Trebitsch et al. 2021). As OBELISK is a zoom simulation, we only considered halos that do not contain any low-resolution DM particle to avoid artefacts. In simulated high-redshift galaxies, disturbed morphologies are common, which makes it challenging to define the centre of a galaxy. We followed what has been done for the NEW-HORIZON simulation (Dubois et al. 2021) and chose as our fiducial ‘centre’ the position of the density peak (for stars, DM, and both) determined recursively using a shrinking sphere approach.

Since BHs are not artificially pinned to galaxy centres, we have to assign BHs to galaxies. The main BH of a galaxy was defined as the most massive BH located within max(4Δx, r₅₀), where r₅₀ is the half-mass radius r₅₀. BHs that are not assigned to any galaxy as main BHs can be assigned as satellite BHs to the highest stellar mass galaxy enclosing them within 2r₅₀. Finally, star formation rates were averaged over 10 Myr. We note that galaxy properties were stored in snapshots, recorded every 15 Myr or less. BH properties were instead recorded at every coarse timestep, about 0.1 Myr.

2.3. Selection of BH mergers

BHs that are merged in the simulation following the sub-grid algorithm for BH-BH mergers (at a distance of 4Δx) were identified as ‘numerical mergers’. To find the mergers host galaxy, we identified in the snapshot immediately after the merger which galaxy contains the location of the merger within a distance < r₅₀ from the galaxy centre. BH mergers occurring at larger distances from all galaxy centres were considered spurious cases and discarded. To account for the continued dynamical decay of the BH binary below 4Δx we calculated delays in post-processing, and defined ‘delayed mergers’ as the outcome of adding such delays. For numerical mergers, the BH properties are measured at the coarse timestep immediately preceding, and the galaxy properties at the simulation snapshot prior to the merger. For numerical and delayed remnants, we measured galaxy properties at the first available post-merger output and BH properties at the first available post-merger coarse timestep.

Sub-grid merger delays were modelled as in Volonteri et al. (2020) and Dong-Páez et al. (2023). We included a dynamical friction phase from the position where the BHs were located at the numerical merger down to the centre of the host galaxy by computing the dynamical friction timescale for a massive object in an isothermal sphere, considering only the stellar component of the galaxy and including a factor 0.3 to account for typical orbits being non-circular. We calculated the sinking time of both BHs in the numerical merger and took the longest of the two. For binaries whose dynamical friction timescale ends before the simulation is stopped at z = 3.5 we further calculated binary evolution timescales through interaction with stars Sesana & Khan (2015) and gas Dotti et al. (2015) until gravitational waves take over. Delayed mergers predicted to occur at z < 3.5, after the final redshift down to which OBELISK has been run, cannot be modelled as we lack the information on host galaxy properties required to compute the relevant timescales.

We note that during this post-processed dynamical evolution, the BHs have already been merged numerically in the simulation. That is, the simulation does not track the individual evolution of the two BHs during the merger delays, but only that of a numerically merged BH with the total mass. Consequently, for delayed mergers, merger parameters that require individual properties of the BHs, such as the mass ratio or the pre-merger BH spins, cannot be extracted from the simulation. These parameters would need to be estimated in post-processing, and the final value would be strongly dependent on the model used (e.g., Farris et al. 2014; Duffell et al. 2020; Muñoz et al. 2020; Siwek et al. 2020). Since mass ratio is a key parameter for GW analysis, we excluded delayed mergers from any analysis involving the mass ratios or spins at merger. We considered both numerical mergers and delayed mergers as a way to bracket our uncertainty.

2.4. AGN spectral energy distribution

Commonly, the AGN luminosity in different bands is estimated using bolometric corrections (e.g., Shen et al. 2020), which are derived from mean observed quasar Spectral Energy Distributions (SEDs). Due to observational selection effects, the quasar samples used to calibrate such models tend to cover only a reduced region of the BH parameter space. In contrast, our simulated BH sample spans a much wider range in M_• and f_Edd.

In order to capture qualitatively the effect of BH physical parameters for a large region of the parameter space, it is often preferable to model AGN emission by adopting physically motivated analytical models (e.g., Kubota & Done 2018). These models should converge to the standard SEDs used to calculate bolometric corrections when they are restricted to the range of M_• and f_Edd that characterise the quasar samples used to calibrate the bolometric corrections (see for instance the Appendix in Volonteri et al. 2017). We modelled the AGN SED as the sum of the emission from a self-gravitating relativistic, geometrically thin, optically thin accretion disc (Novikov & Thorne 1973) and a power-law X-ray emission with an exponential cutoff from the corona,

$$\begin{array}{c}\hfill L_{\nu }=\{\begin{array}{cc}\mathcal{A}{L}_{\nu ,\mathrm{disc}}\hfill & \nu \le {\nu }_{\max }\hfill \\ \mathcal{A}{L}_{\nu ,\mathrm{disc}}+\mathcal{B}{L}_{\nu ,\mathrm{disc}}(c/2500\mathrm{\AA }){\left(\frac{\nu }{(c/2500\mathrm{\AA })}\right)}^{{\alpha }_{\mathrm{X}}}{e}^{-\nu /{\nu }_{\mathrm{c}}}\hfill & \nu >{\nu }_{\max }\hfill \end{array}\end{array}$$(3)

where 𝒜 and ℬ are normalisation constants and L_ν, disc is the emission from a Novikov-Thorne accretion disc. In Eq. (3), the first term corresponds to the thermal disc contribution, and the second term corresponds to a power-law emission dominating at high energy. The second term is switched on at ν_max the frequency at which the Novikov-Thorne solution peaks. We fix the power-law index of the second term to α_X = −0.9 and the characteristic cutoff frequency to ν_c = 300 keV/h (e.g., Shen et al. 2020).

The disc is assumed to radiate as a blackbody at each radius, with a temperature given by the Novikov-Thorne solution (Krolik 1999). The total disc SED was then obtained by integrating the blackbody emission of each annulus over the radial extent of the disc. The disc is assumed to extend from the radius of the innermost stable circular orbit r_isco. We set the maximum radius of the disc where the disc becomes self-gravitating. For a radiation pressure-dominated disc with opacity given by electron scattering, this radius is given by (Laor & Netzer 1989)

$$\begin{array}{c}\hfill r_{\mathrm{sg}}\simeq 2150{\alpha }^{2/9}{\left(\frac{M_{\bullet }}{{10}^9{M}_{\odot }}\right)}^{-2/9}{\dot{m}}^{4/9}.\end{array}$$(4)

The viscosity parameter α was set to 0.1 (in agreement with numerical studies, e.g., Hawley & Krolik 2002; Hirose et al. 2009) and ṁ ≡ (ε_r(0)/ε_r(a)(L/L_Edd).

The normalisation constants 𝒜 and ℬ were set so that (i) the total integrated luminosity equates to the AGN bolometric luminosity L = ε_r(a)Ṁ_∙c² and (ii) the relative normalisation between the optical and X-ray luminosities, which is characterised by the parameter α_OX,

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{OX}}=0.384{log}_{10}\left(\frac{L_{\nu }(2\mathrm{keV}/h)}{L_{\nu }(c/2500\mathrm{\AA })}\right),\end{array}$$(5)

fits the physical models by Done et al. (2012) and Dong et al. (2012), which roughly predict

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{OX}}=-0.13{log}_{10}(M_{\bullet }/M_{\odot })+0.15{log}_{10}(f_{\mathrm{Edd}})-0.45.\end{array}$$(6)

In practice, the spectrum at 2500 Å and 2 keV is dominated respectively by the disc and the corona, which implies that the ratio between the normalisation constants is approximately given by $\mathcal{A}/\mathcal{B}\approx {10}^{{\alpha }_{\mathrm{OX}}/0.384}{(2500\mathrm{\AA }/2\mathrm{keV})}^{{\alpha }_{\mathrm{X}}}$. Finally, the values of 𝒜 and ℬ can be obtained by integrating Eq. (3) over frequency and imposing condition (i) above. We did not consider a reduction in the flux from the disc due to the random viewing angle. However, the correction would be on average of order unity and would only affect the UV fluxes.

In summary, this model depends on three BH parameters, M_•, f_Edd and a, with only a weak dependence on the spin. A simplified version of this model was used in Trebitsch et al. (2021) to calculate the ionising radiation in OBELISK, therefore ensuring consistency between the in-simulation AGN properties and those calculated in post-processing. This SED is generally appropriate for radiatively efficient discs, which characterise BHs with f_Edd ≳ f_Edd, crit = 0.01. Otherwise, it can be regarded as an upper limit. We have checked that in our sample all BHs detectable in UV or X-ray and more than 94% of those detectable in radio have f_Edd > f_Edd, crit (see Sects. 3.2 and 3.3).

We used this model to calculate the flux density in the rest-frame UV at λ = 1500 Å, which corresponds to optical to near-IR observed emission taking the redshift of sources, z > 3.5, into account. This is therefore the observability with optical telescopes. We calculated the integrated flux in the observer-frame ‘total’ X-ray band [0.5 − 10] keV. We further define the observer-frame soft, [0.5 − 2] keV, and hard, [2 − 10] keV, X-ray bands.

From a differential luminosity L_ν, the spectral flux density at an observed frequency ν can be calculated from luminosity at the rest-frame frequency ν_em = (1 + z)ν,

$$\begin{array}{c}\hfill S_{\nu }=(1+z)\frac{L_{{\nu }_{\mathrm{em}}}}{4\pi D_L^2(z)}\end{array}$$(7)

where the (1 + z) factor reflects the redshifting of the differential bandwidth dν. The cosmological luminosity distance D_L(z) of the source was calculated using the fiducial cosmology of the OBELISK simulation (see Sect. 2.1).

The integrated flux in a given band can be calculated from the integrated rest-frame luminosity L as follows,

$$\begin{array}{c}\hfill F=\frac{L}{4\pi D_L^2(z)}.\end{array}$$(8)

We show an example of this SED model in Fig. 2 for a merger remnant with M_• = 1.3 × 10⁸ M_⊙, f_Edd = 0.51 and a = 0.90, at z = 3.55, and compare it with the emission from the host galaxy and some realistic instrumental sensitivity limits, which are defined in the sections below.

Fig. 2. Spectral energy distribution of a merger remnant with mass 1.3 × 108 M⊙, fEdd = 0.51 and spin 0.90, residing in a galaxy with mass 1.4 × 1011 M⊙ at redshift 3.55. The top panel depicts the (unobscured) rest-frame luminosities while the bottom panel depicts the (obscured) observer-frame flux. The SED of the remnant BH is shown in the solid black and red lines for the fiducial and the brightening scenario, while the approximate UV and X-ray galaxy SED is shown in the blue dashed lines. The green arrows indicate approximately the instrumental sensitivities assumed in radio at 2 GHz (assuming a sensitivity of 0.4 μjy), in UV at 1500 Å (mUV, lim = 31.3), and in the 0.5 − 2 keV band (FX, lim = 10−17 erg s−1).

2.5. AGN obscuration

We estimated the column density of gas in the interstellar medium (ISM) contributing to the BH obscuration by casting rays isotropically around each BH in the outputs of the simulation. For each BH we casted 100 rays and integrated the gas column density from the BH position to the virial radius of its host halo. We used a version of the RASCAS code (Michel-Dansac et al. 2020) modified to integrate column densities very efficiently. For each sightline i, we computed exp(−N_H, iσ_T) and defined the typical column density around each BH as

$$\begin{array}{c}\hfill {\overline{N}}_{\mathrm{H}}=-ln\left(\langle exp(-N_{\mathrm{H},i}{\sigma }_{\mathrm{T}})\rangle \right).\end{array}$$(9)

Additionally, we incorporated a subgrid model in order to account for the unresolved contribution from the pc-scale, geometrically thick gas structure surrounding accreting BHs generally known as the torus. We assumed that the gas is at rest at infinity and that it falls radially under the action of the BH’s gravitational pull. This assumption gives a lower limit on the time spent inside the torus by the accreted material, and therefore a lower limit on N_H. We assumed a spherically symmetric configuration, which would roughly correspond to a spherically averaged problem. Under these assumptions, the gas density in the torus should follow a power-law density profile ρ ∝ r^−3/2. We normalised the density so that the total mass in the torus, between its inner (r_in) and outer (r_out) radii is equal to Ṁ_∙t_ff. This is because the material crosses the torus in a free-fall time (Hönig & Beckert 2007) $t_{\text{ff}}\sim \pi {(r_{\text{out}}^3/(8GM))}^{1/2}$, where we have assumed r_in ≪ r_out. The inner radius can be assumed to be the dust sublimation radius (Suganuma et al. 2006),

$$\begin{array}{c}\hfill r_{\mathrm{in}}\approx 0.47{\left(\frac{L_{\mathrm{UV}}}{{10}^{46}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{1/2}\mathrm{pc}.\end{array}$$(10)

The outer radius can be approximated as the BH gravitational sphere of influence for gas, that is, the BHL radius, at $r_{\mathrm{out}}=2GM/({\overline{c_{\mathrm{s}}}}^2+{\overline{v}}_{\mathrm{rel}}^2)$.

Further assuming that the infalling gas is only composed of hydrogen, we integrated the density in the radial direction to arrive at the following expression

$$\begin{array}{c}\hfill {\overline{N}}_{\mathrm{H},\mathrm{torus}}=\frac{3\pi }{4{\sigma }_{\mathrm{T}}}\frac{f_{\mathrm{Edd}}}{{\epsilon }_{\mathrm{r}}}\frac{{\overline{c_{\mathrm{s}}}}^2+{\overline{v}}_{\mathrm{rel}}^2}{c}(\sqrt{\frac{r_{\mathrm{out}}}{r_{\mathrm{in}}}}-1).\end{array}$$(11)

We note that the equation depends mainly on the Eddington ratio f_Edd, but depends also on the BH mass, the radiative efficiency ε_r, and the local gas properties.

In the X-rays, the rest-frame attenuated luminosity can be calculated as

$$\begin{array}{c}\hfill L_{\mathrm{X},\mathrm{abs}}={\displaystyle {\int }_{{\nu }_{\min }}^{{\nu }_{\max }}}{L}_{\nu }{e}^{-{\sigma }_{\mathrm{X}}(\nu )N_{\mathrm{H}}}\mathrm{d}\nu .\end{array}$$(12)

The X-ray cross sections are calculated from the polynomial fits in Morrison & McCammon (1983), extrapolated if needed to ν > 10 keV h⁻¹ assuming a scaling σ_X ∝ ν⁻³.

Given that OBELISK includes a model for dust, we followed a similar approach to estimate the UV obscuration. For each BH, we casted 100 rays in different directions. Along each sightline i, we computed the dust optical depth τ_UV, i = κ₁₅₀₀∫_lρ_d, idl where ρ_d, i is the local mass density of dust in each cell along the sightline i and κ₁₅₀₀ is the dust mass absorption coefficient at 1500 Å. We estimated κ₁₅₀₀ as follows: we started by assuming that our dust is composed of a mixture of silicate and carbonaceous grains with respective mass fractions 54% and 46% as in Aoyama et al. (2018), following Hirashita & Yan (2009), and that the grain size distribution follows the MRN grain size distribution (Mathis et al. 1977) between 0.001 and 0.25 μm. We then integrated the extinction cross sections from Laor & Draine (1993) over the grain size distribution to get κ₁₅₀₀. The typical attenuation is then defined as previously by

$$\begin{array}{c}\hfill {\overline{\tau }}_{\mathrm{UV}}=-ln\left(\langle exp(-{\tau }_{\mathrm{UV},i})\rangle \right).\end{array}$$(13)

We similarly added a torus correction to the UV obscuration. In order to obtain a torus UV optical depth, we multiplied N_H, torus by a factor κ₁₅₀₀M_dust, 50/M_gas, 50, where M_dust, 50 and M_gas, 50 are the dust and gas mass inside r₅₀.

The median value and interquartile scatter of the total gas column density for our numerical merger population is log₁₀(N_H/cm⁻²) = 23.3 ± 0.2. The interstellar contribution (log₁₀(N_H, ISM/cm⁻²) = 23.3 ± 0.3) dominates, while the torus contribution (log₁₀(N_H, torus/cm⁻²) = 22.2 ± 0.4) only accounts for less than 10% of the total median value. For dust, the median UV optical depth is log₁₀τ_UV = 0.97 ± 0.28. In this case, the interstellar contribution still dominates (log₁₀τ_UV, ISM = 0.84 ± 0.31), but the torus contributes more significantly (log₁₀τ_UV, torus = 0.11 ± 0.46), to 30% of the median value. We do not consider the contribution from the intergalactic medium since it is subdominant for our sample (Arcodia et al. 2018).

We can define the mean optical depths in the observer-frame soft and hard X-rays τ_Xs and τ_Xh as τ = ln(F/F_abs), where F and F_abs are the integrated unabsorbed and absorbed flux. The distributions of τ_UV, τ_Xs, and τ_Xh as a function of host galaxy mass M_* are shown in Fig. 3 for remnant BHs in our sample of numerical mergers. The optical depths are generally very high in the UV. In the X-rays, the obscuration is much smaller, since at the high rest-frame frequencies probed the gas cross-sections are very small. Recall the simulation is limited at high redshift, so all our mergers occur at z ≳ 3.5.

Fig. 3. Distribution of the optical depth against the host galaxy stellar mass for remnant BHs in our numerical merger sample. Top panel: rest-frame UV; middle panel: observer-frame 0.5 − 2 keV; bottom panel: observer-frame 2 − 10 keV. The colour denotes the number of BHs in each bin, on a logarithmic scale. Optical depths are very high in the UV, while in X-rays the optical depths are not extreme.

The effect of obscuration can be seen in the bottom panel can be seen in Fig. 2, where we show an example of an obscured observer-frame SED. The effect of obscuration is particularly strong in the UV.

2.6. Radio emission

In order to assign a jet radio luminosity to the simulated mergers while bypassing the theoretical uncertainties related to the production of jets, we resorted to the ‘fundamental plane of black hole activity’, an empirical correlation between the radio luminosity, X-ray luminosity, and mass of BHs. This relation has been shown to be applicable to BHs spanning 8 orders of magnitude in mass (e.g., Merloni et al. 2003; Falcke et al. 2004). Gültekin et al. (2014) proposed that the fundamental plane relation found in Gültekin et al. (2009) also holds for a sample of low-mass (M_• ≲ 10^6.3 M_⊙), highly accreting AGNs. More recent work (Gültekin et al. 2022) has however shown that the fundamental plane tends to underestimate L_R at fixed L_X for a sample of highly accreting AGN powered by low-mass BHs. The fundamental plane only takes into account the contribution of the core radio luminosity. Therefore, we computed radio luminosity from the relation found in Gültekin et al. (2009), but we treat this ‘pessimistic’ model as a lower limit in the following analysis. We calculated the radio luminosity as

$$\begin{array}{c}\hfill {log}_{10}{L}_{\mathrm{R},5\mathrm{GHz}}=4.80+0.78{log}_{10}{M}_{\bullet }+0.67{log}_{10}{L}_{\mathrm{X}},\end{array}$$(14)

where L_R, 5 GHz = ν_5 GHzL_5 GHz is the radio luminosity at 5 GHz, and L_X is the integrated X-ray luminosity in the rest-frame 2 keV − 10 keV energy range. We calculated the radio luminosity at 2 GHz assuming a power-law spectrum L_ν ∝ ν^α_R with index α_R = −0.7 (Gültekin et al. 2014).

As an upper limit to the radio luminosity, we considered a theoretical model in which the jet is powered by the Blandford-Znajek effect (Blandford & Znajek 1977). We modelled the total synchrotron luminosity based on Meier (2001),

$$\begin{array}{c}\hfill L_{\mathrm{S},\mathrm{s}}=\{\begin{array}{cc}\begin{array}{c}{10}^{46.0}{\eta }_{\mathrm{S}}{\left(\frac{{\alpha }^{\mathrm{AD}}}{0.3}\right)}^{-1}\left(\frac{M_{\bullet }}{{10}^9{M}_{\odot }}\right)\times \hfill \\ \left(\frac{f_{\mathrm{Edd}}}{0.1}\right)g^2(0.55f^2+1.5fa+a^2)\mathrm{erg}{\mathrm{s}}^{-1}\hfill \end{array}\hfill & f_{\mathrm{Edd}}<0.01\hfill \\ \begin{array}{c}{10}^{43.5}{\eta }_{\mathrm{S}}{\left(\frac{\alpha }{0.01}\right)}^{-0.1}{\left(\frac{M_{\bullet }}{{10}^9{M}_{\odot }}\right)}^{0.9}\times \hfill \\ {\left(\frac{f_{\mathrm{Edd}}}{0.1}\right)}^{1.2}(1+1.1a+0.29a^2)\mathrm{erg}{\mathrm{s}}^{-1}\hfill \end{array}\hfill & f_{\mathrm{Edd}}\ge 0.01.\hfill \end{array}\end{array}$$(15)

Here, the top equation represents the jet production from a geometrically thick advection-dominated accretion flow (ADAF) at low f_Edd, while the bottom equation corresponds to the thin disc case. As above, we assumed for the thin disc α = 0.1, while for the ADAF we assumed α = 0.3 following Meier (2001). We set the dimensionless constant f and g to 1 and 2.3, following Meier (2001) and Tamanini et al. (2016). We further assumed, following Meier (2001) that only a fraction η_S = 10⁻² of this power is radiated in the synchrotron spectrum. In contrast to the fundamental plane, this expression provides the total jet power, not only the core jet power.

The radio luminosity νL_ν at ν = 2 GHz can be calculated assuming the synchrotron radiation is emitted in a power-law spectrum with index α_R over the frequency range 0.01 − 100 GHz. Overall, a fraction ∼10⁻³ of the initial power is transformed into radio luminosity νL_ν at 2 GHz. In general, Eq. (15) can predict radio luminosities more than 2 orders of magnitude above the fundamental plane. Thus, we regard this ‘optimistic’ model as an upper limit in the following analysis.

We show an example of the optimistic radio model in Fig. 2.

2.7. Merger-induced brightenings and transients

A BH merger can potentially induce a brightening on a scale of days to years around the time of merger, which increases the luminosity of the remnant BH and constitutes a transient signal that can be used to detect the merger and identify it as such. To model this, we assumed that initially the remnant BHs emits in X-ray and UV at the fiducial luminosity predicted by our SED model above, for the appropriate accretion rate calculated in the simulation. That is, we assumed that the SED model for a single BH in an α disc applies. A brightening occurs either shortly before the merger (Armitage & Natarajan 2002; Cerioli et al. 2016) or after ∼t_vis, when the accretion of the inner rim drives a burst. We remain agnostic on the exact process, and we modelled a brightening by assuming f_Edd = 1 in our SED model¹. Physical arguments suggest that the formation of the cavity and the subsequent brightening should only happen if q is large enough and if the binary is embedded in a gas-rich environment, but we optimistically assume this to apply to all mergers. On the other hand, we did not include a pre-merger suppression of the accretion rate, which is a pessimistic approach, in the sense that the luminosity change is weaker than if accretion were suppressed.

Different changes to jet properties around the time of BH mergers have also been proposed. Here, we considered the possibility of an increase in the jet radio luminosity due to an increase in the accretion rate analogous to the UV and X-ray model above. Again, we modelled this brightening by assuming f_Edd = 1 in our radio models.

We also considered a short-lived flare as a transient feature and modelled it as having, for an equal mass merger, a luminosity a factor (1 + 𝒦) higher than the pre-merger luminosity L_S, s1 + L_S, s2. The increase in luminosity due to a flare, parametrised here by 𝒦, is not well constrained by simulations (Moesta et al. 2012; Gold et al. 2014; Kelly et al. 2017; Cattorini et al. 2021, 2022), and can depend strongly on the merger parameters. We assumed a factor of 𝒦 = 5. Further, following the discussion in Kaplan et al. (2011), we assumed an approximate scaling of the flare luminosity $L_{\text{S},\text{f}}\propto q^2{M}_{•}^2$. For simplicity, we neglected the dependence on other parameters. In summary, we assumed that

$$\begin{array}{c}\hfill L_{\mathrm{S},\mathrm{f}}=(1+\mathcal{K}{q}^2)(L_{\mathrm{S},\mathrm{s}1}+L_{\mathrm{S},\mathrm{s}2}).\end{array}$$(16)

These two merger-induced transient phenomena (brightening and flare) will alter the small-scale core radio luminosity, but not necesarily the total radio luminosity. Therefore, we apply this radio transient model only on the ‘pessimistic model’, which estimates only the core radio emission.

In our model, we did not consider the possibility explored by Yuan et al. (2021) and Ravi (2018) that a gamma-ray burst-like source is produced as a newly formed jet after the merger impacts and shocks the nearby gas. This choice is based on simulations showing that the jet can exist both before and after the merger. Even in the case of a spin flip, Kelly et al. (2021) find that the jet direction does not change, remaining aligned with the ambient magnetic field on large enough scales. Similarly, Ruiz et al. (2023) do not find a significant perturbation in the jet propagation due to the change in spin direction for initially slowly spinning BHs. This presumably means that the BH will eventually realign with the jet, although the outcome is unclear and depends on complex physics (see McKinney et al. 2013; Liska et al. 2021). In this case the jet will continue to propagate in the same direction as before, rather than encountering pristine gas.

Finally, we note that other types of transient features can appear around the time of the merger: spectral changes caused by perturbations in the accretion disc (Schnittman & Krolik 2008), complex lightcurves in the case of kicked BHs (Rossi et al. 2010), periodic modulations (Gutiérrez et al. 2022). Exploration of these and also of features occurring during the inspiral (Sesana et al. 2012; Roedig et al. 2014; Farris et al. 2015) are postponed to a future investigation.

2.8. Galactic emission

We also modelled the galactic emission, which in our case acts as contamination hindering the detection of the BH merger. In the UV band, the galactic emission was computed from the stellar population in the galaxy. For each galaxy in our catalogue, we estimated the intrinsic UV luminosity at 1500 Å from the properties of the star particles in the galaxy. Each star particle was attributed a luminosity based on its age and metallicity using the BPASS v2.2.1 SED (Eldridge et al. 2017; Stanway & Eldridge 2018) and rescaled to the mass of the star particle. The intrinsic luminosity of the galaxy was then obtained by summing over all star particles associated with the galaxy. As OBELISK includes a model for the formation, growth, and evolution of dust in the galaxy, we can estimate the observed UV luminosity by computing the average attenuation. For this, we used the RASCAS code to cast 100 rays isotropically from each star particle within 5r₅₀ of each galaxy. We measured the dust optical depth along each ray and used the average f_UV = ⟨exp(−τ_UV, i)⟩ as the escape fraction of UV light. The observed UV luminosity of a galaxy was then defined as the intrinsic luminosity times f_UV. We note that our method does not account for orientation effects.

The galactic X-ray emission was assumed to be dominated by X-ray binaries (XRBs). We modelled the X-ray luminosity of XRBs using the empirical scaling relation found in Fornasini et al. (2018), in which the integrated X-ray luminosity in the 2 keV–10 keV range is parametrised as a function of the galaxy stellar mass M_* and star formation rate (SFR) as follows,

$$\begin{array}{cc}\hfill L_{\mathrm{XRB}}& ={10}^{29.98}{(1+z)}^{0.62}\left(\frac{M_{\ast }}{M_{\odot }}\right)\mathrm{erg}{\mathrm{s}}^{-1}\hfill \\ \hfill & +{10}^{39.78}{(1+z)}^{0.2}{\left(\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}\right)}^{0.84}\mathrm{erg}{\mathrm{s}}^{-1}.\hfill \end{array}$$(17)

The effect of gas absorption was added assuming a constant column density of N_H = 10²² cm⁻² and a power-law spectrum with photon index Γ = 1.4.

We compare the host galactic emission with the remnant BH of a particular merger in Fig. 2. In general, the galactic emission can be comparable with the AGN emission.

The radio emission generated by star-forming regions can also hinder the radio detectability of AGN. As an order of magnitude estimate, we calculated the galactic radio emission from the scaling relations in Bell (2003), which relate the radio luminosity L_R, gal to the SFR. Since those estimates are based on a fit of SFR as a function of L_R, gal and not the converse, which is what we need, and the data is limited to low redshifts, we do not include them explicitly in our analysis, but we note that the SFR-induced radio luminosities could be comparable or higher than the AGN for BHs with mass M_• ≲ 10^7.5 M_⊙.

2.9. LISA gravitational wave analysis

We calculated the detectability and parameter estimation for the satellite LISA for our set of simulated mergers. Since delayed mergers do not have a well-defined mass ratio, as it is unclear how much mass is accreted onto which black hole during the sub-grid inspiral, we restricted our GW analysis to the sample of numerical mergers.

To simulate the GW signal, we adopted the PhenomHM waveform (London et al. 2018) that assumes spins aligned with the binary orbital momentum but includes higher order harmonics to break degeneracies in the parameter estimation process. The GW signal from a BH-BH binary with aligned spins can be described by 11 parameters: the primary and secondary mass M₁ and M₂ (with M₁ > M₂), the component of the spins aligned to the orbital angular momentum χ₁ and χ₂ (${\chi }_{1,2}={\mathit{a}}_{\mathbf{1},\mathbf{2}}.\mathit{\ell }/\ell$), the time of coalescence t_c, the luminosity distance D_L(z), the inclination ι, the sky latitude β and longitude λ, the orbital phase at coalescence ϕ, and the polarisation ψ. We can additionally define the chirp mass M_chirp = (M₁M₂)^0.6(M₁ + M₂)^−0.2. M₁, M₂, χ₁ and χ₂ are directly produced in the simulation. The sky latitude β and longitude λ were set randomly over the sphere as well as the inclination ι and the polarisation ψ. The phase at coalescence ϕ is randomised between [0, 2π]. The time to coalescence t_c was set randomly between [0, 1] yr.

For a single event, we computed the signal-to-noise ratio (S/N) as

$$\begin{array}{c}\hfill S/N={\displaystyle {\int }_{{\nu }_{\mathrm{lo}}}^{{\nu }_{\mathrm{up}}}\frac{|\stackrel{\sim }{h}{(\nu )|}^2}{S_n(\nu )}\mathrm{d}\nu ,}\end{array}$$(18)

where $\stackrel{\sim }{h}(\nu )$ is the Fourier transform of the time-domain GW signal, S_n(ν) is the noise power spectral density, and ν_lo and ν_up are the minimum and maximum frequencies of integration. For S_n(ν) we adopted the so-called ‘SciRDv1’ sensitivity (Babak et al. 2021) and we set ν_lo = 10⁻⁵ Hz and ν_up = 0.5 Hz. If t_c is very short, the initial frequency ν_lo was reduced. We also added to the LISA power spectral density the noise from the population of unresolved galactic binaries (Karnesis et al. 2021), with an amplitude corresponding to three years of observations. The posterior distributions on the binary parameters $\overline{\theta }$ can be obtained following Bayes theorem as

$$\begin{array}{c}\hfill p(\overline{\theta }|d)=\frac{\mathcal{L}(d|\overline{\theta })\pi (\overline{\theta })}{p(d)},\end{array}$$(19)

where $\mathcal{L}(d|\overline{\theta })$ is the likelihood of observation d with parameters θ, $\pi (\overline{\theta })$ are the prior probabilities on the binary parameters and $p(d)=\int \mathrm{d}\overline{\theta }\mathcal{L}(d|\overline{\theta })\pi (\overline{\theta })$ is the evidence.

The binary posterior distributions were obtained following the formalism presented in Marsat et al. (2021). For each binary, we ran the Bayesian Markov chain Monte-Carlo (MCMC) analysis for 2000 iterations with 64 walkers and 10 temperatures. We added an additional step to our analysis. Since we are interested in the sky localisation Δ#x03A9; to detect the possible EM emission, we decided to re-run the systems with 5 < Δ#x03A9;/deg² < 40 for 10⁵ iterations to ensure the convergence of the MCMC algorithm.

3. Multi-messenger observability of BH mergers

In this section, we analyse the detectability of merging BHs at all redshifts probed by the simulation, z ≥ 3.5. We analyse the GW emission of the sample of numerical mergers. We recall that we consider only numerical mergers because the mass ratio q, which is crucial for GW data analysis, is not well-defined for delayed mergers (see Sect. 2.3). The analysis of EM signals is extended to both numerical and delayed mergers.

We define a number of samples that are favourable for an EM detection. We define a source as AGN-dominated if the flux from the BH is larger than the flux from its host galaxy, and so the galactic contamination does not hinder the detection of the BH. We denote the converse case, where the galactic flux dominates the BH flux, as galaxy-dominated. A BH is considered observable if its flux exceeds the instrument sensitivity and the source is AGN-dominated. If the GW analysis is performed for the sample, we additionally require that a GW signal be detectable by LISA to consider it observable as a multi-messenger source. We note that this definition does not require an EM transient to be present.

We define an EM counterpart to a GW event as a source that exhibits a merger-induced transient with a significant change in flux, enabling identification as a merger. We assume that such a merger-induced variation in the flux can be detected if either: (i) The source ‘appears’ at the time of the transient – it is undetected before the transient and detected at the transient. The BH is observable at the transient. (ii) The source ‘disappears’ – it is detected before and undetected after. The BH is observable before. (iii) The source is detected before and after, but the flux changes significantly. In this case (iii), we consider a flux difference to be significant enough if the flux varies by more than a factor of 2. This is probably an optimistic choice since AGN are intrinsically variable sources. Additionally, the BH must be observable either before or at the transient. We denote mergers that fall into any of these categories as having an EM counterpart. In the following, we do not consider transients of type (ii) since in the analysis below we do not find any transients of this type.

3.1. GW observability and LISA parameter estimation

We first discuss the detectability and parameter estimation by the LISA satellite. The analysis is performed by accumulating signal from the time the binary enters the LISA band up until coalescence. Figure 4 shows for the numerical merger population the distribution of the S/N, and the distribution of 90%-confidence uncertainties in estimating z, M_chirp, q, χ₁ and χ₂, and the sky localisation. We find that most of the mergers in our sample are detectable if we set the limit for detectability at S/N > 10. This is because the masses of our BH population generally lie within the target range of LISA, and the mass ratios are generally moderate. An exception is the lowest mass ratio mergers (q ≲ 10⁻³), which in our simulation correspond to mergers with a total mass of M_• ≳ 10⁸ M_⊙. These mergers have low S/N despite occurring at the lowest redshifts of our sample, z < 3.7. Additionally, some high-redshift (z ≳ 4) edge-on mergers of BHs close to the seed mass (M_chirp ≲ 10⁵ M_⊙) are undetected or detected with poor S/N.

Fig. 4. Parameter distribution for the mock LISA data analysis. The top-left panel shows the signal-to-noise ratio (S/N) distributions. Gravitational wave events are assumed to be detected if S/N > 10, indicated by the vertical dashed line. The 90%-confidence error in the sky localisation is shown in the bottom-left panel. The dashed line delimits the shaded region, which corresponds to the events with sky localisation poorer than 10 deg2. The relative 90%-confidence error distributions for z, Mchirp and q, and for χ1 and χ2, are shown in the right panels, top and bottom respectively. Most parameters are recovered with small uncertainties. The majority of mergers have sky localisation worse than 10 deg2 at the time of the merger.

The parameters encoding the binary masses, M_chirp and q are estimated with high precision, especially M_chirp. The redshift and the spins are recovered with good precision, but the distributions have tails extending to high uncertainties that may become comparable to the value of the parameters. The 90%-confidence uncertainties in the sky localisation at merger, which corresponds to the best available estimate (compared to considering the inspiral phase only, see Mangiagli et al. 2020), are generally poor – only 37% of the mergers have a localisation better than 10 deg². The redshift determination has typical uncertainty of 0.01–0.1, which means that the 3D error box is dominated by the sky localisation uncertainty. As a reference, in COSMOS2020 (Weaver et al. 2022) at i-magnitude < 27 and 3.49 < z < 3.51 there are 740 galaxies per deg² (Shuntov, priv. comm.).

This hinders the possibility of using GW detection to guide the search for EM counterparts with instruments having a small field of view. In radio, the field of view of ngVLA and SKA are below 2 deg² so in most cases several pointings are needed in order to cover the LISA sky localisation error-box. In X-ray the field of view of Athena is planned to be 0.4 deg² (0.25 deg² for NewAthena) and that of AXIS has a proposed ∼0.13 deg², which will require more tiling, while the NASA Transient Astrophysics Probe is proposed to be 1 deg², but with a lower sensitivity than Athena and AXIS. THESEUS has a very large field of view, 0.5–2 sr, but it is expected to have much lower sensitivity. In optical, with large field-of-view instruments such as the Rubin Observatory, 9.6 deg² (Ivezić et al. 2019), one can use only a few tiles to cover the error-box.

The uncertainty in the sky localisation of our systems is mainly determined by the inclination ı, the angle of the orbital angular momentum of the binary with respect to the line of sight. Face-on systems are better localised, leading to an error of ∼10⁻¹ on average, while edge-on systems lead to ∼10³. Since our systems are distributed uniformly in orientation, the inclination angle ı is randomly distributed with a probability proportional to sin ı. This results in a large scatter in the sky localisation with values preferentially skewed towards the poorly localised regime. It is important to note that the waveforms used in our parameter estimation do not include the effects of spin precession. The inclusion of spin precession could improve the errors by a factor of 2 − 5 (Lang & Hughes 2006), which in our case would lead to 46 − 60% of the mergers being localised better than 10 deg² if all errors were scaled down equally. Our current value of 37% is likely a lower limit.

If merger parameters are estimated from the median of our parameter distributions, we generally recover the true values with good accuracy. However, it is worth noting that there is a small number of extreme outliers in our sample. These are generally events detected with low S/N, for which the parameter estimation returns large errors and uncertainties. These outliers tend to strongly overestimate the luminosity distance (reaching values of up to z ∼ 40 assuming the fiducial cosmology) and chirp mass of the events, which may be interpreted erroneously as evidence for massive primordial black holes (see De Luca et al. 2021; Ng et al. 2022; Martinelli et al. 2022, for detailed models of how GWs can constrain primordial black holes). In all these cases, the high deviations are accompanied by higher uncertainties, which offer a way to flag them as outliers.

3.2. UV and X-rays

3.2.1. UV and X-ray detectability

In the following sections, we explore the possibility of an EM detection that would complement the GW detections discussed in the previous section. In Fig. 5, we show for numerical and delayed merger remnants the attenuated remnant AGN flux in the rest-frame UV (1500 Å) and the observer-frame 0.5 − 10 keV band X-rays against the host galaxy flux, which in our case acts as a contaminant hindering the detection of the central AGN. The rest-frame UV wavelength considered corresponds to the optical or near-infrared in the observer frame for our high-redshift sample – at z = 3.5 it would be observed in the g-band, while at z = 7 it would be observed in the z-band. For the modelling of the AGN SED, the implementation of obscuration, and the estimation of the galactic emission, we refer to Sects. 2.4, 2.5, and 2.8, respectively. For our merger hosts, the XRB luminosity is generally dominated by high-mass XRBs, whose luminosity depends on the SFR.

Fig. 5. Black-hole vs galaxy flux, both corrected for absorption. Left column: rest-frame UV magnitude at ν = 1500 Å; right column: X-ray 0.5 − 10 keV band. Top row: numerical mergers; bottom row: delayed mergers. Black dots: accretion rate measured in the simulation; red crosses: merger-induced brightening with fEdd = 1. The grey region corresponds to galaxies dominating the emission over their BHs. The black dashed line indicates the boundary where the luminosity and galaxy and BH are equal. The hatched region corresponds to galaxies and BHs below the minimum flux needed for detection, which we set to mUV, lim = 31.3 (UV) and FX, lim = 10−17 erg cm−2 s−1 (X-rays). In the top left corner, we indicate the fraction of BHs with higher flux than their host galaxies and that of BHs above the observational flux limit. We note that a significant fraction of sources are very faint and lie beyond the limits of the plot, especially in the UV. Mergers that are GW undetected by LISA are shown with bigger and lighter markers, in purple and orange, for the fiducial and brightening cases. In UV all merging BHs are outshone by their host galaxy, while in X-rays a fraction of merging BHs are both brighter than the host galaxy and detectable by sensitive X-ray telescopes.

Recall that we refer to sources in which the remnant BH dominates its host as AGN-dominated, while we refer to the converse case where the host dominates over the BH emission as galaxy-dominated. We also set optimistic observability thresholds, at m_UV, lim = 31.3 in the UV (corresponding to the sensitivity for isolated point sources to 5σ in 5 h of the MICADO instrument on ELT², Davies et al. 2010), and F_X, lim = 10⁻¹⁷ erg cm⁻² s⁻¹ in the X-rays, a reasonable flux limit for missions such as AXIS (Mushotzky 2018) and Athena (Nandra et al. 2013). We refer to BHs that exceed this limit and dominate the galactic emission as observable. For numerical mergers, we additionally require the merger to be detectable by LISA for it to be considered observable as a multi-messenger source. We note that our definition of observable does not require a transient to be present.

We find virtually no observable merger remnant BHs in the restframe UV. The few cases that exceed the observability threshold are galaxy-dominated because the high-redshift galaxies in our sample are actively star-forming. Even in the presence of a merger-driven brightening we find that it remains unlikely to detect our sample of merger remnants. If we do not consider dust obscuration, which is very high for AGN in the UV (see Fig. 3), all sources also remain galaxy-dominated.

In X-rays, the detectability of BH mergers is more favourable. A fraction ∼15%−35% (134 out of 878–46 out of 129) of sources are AGN-dominated, depending on whether numerical or delayed mergers are considered. This is because the galactic emission is lower relative to the BH emission in X-rays compared to UV, and gas obscuration is only a small effect due to the small gas cross-sections at the very high energies probed (0.5(1 + z)−10(1 + z) keV in the rest-frame X-ray, see Fig. 3). A fraction ∼5%−14% (42 out of 878–18 out of 129) of all BH mergers is observable. Most BHs above the observability limit also dominate over the galactic emission. Unfortunately, some of the brightest mergers in the X-rays are not observable in GWs. These events correspond to a population of high-mass, low-q mergers, which have poor detectability.

The fraction of observable remnant AGN in X-rays is larger for delayed mergers compared to numerical mergers. Delayed mergers remnants tend to be more massive compared to numerical mergers since the galaxy merger can induce an accretion burst that increases the mass of the BHs and delayed mergers are biased towards high-mass galaxies that can experience mergers early (Dong-Páez et al. 2023). A more massive BH population leads to a higher number of observable mergers. The fraction of AGN-dominated is also higher for delayed mergers, for two reasons: firstly, delayed mergers tend to have higher M_•/M_* ratios than numerical mergers, since delayed merger remnants can be slightly overmassive at fixed M_*. Secondly, the galaxy merger that precedes the BH merger can drive both sSFR and BH accretion rate bursts, which boost the XRB and BH luminosity. In Dong-Páez et al. (2023), we find that the increased sSFR has decayed by the time of the delayed merger, and so the galactic XRB boost is absent for delayed mergers. In contrast, the mass gained by BHs during the BH accretion rate boost allows them to accrete at a higher rate on average even at the delayed merger.

An Eddington-limited brightening can greatly increase the observability of BHs in the X-ray. The fraction of both AGN-dominated and observable sources increases significantly with respect to the fiducial accretion rate, to ∼45%−75% (129 out of 878–98 out of 129) and 10 − 30% (98 out of 878–36 out of 129). We remind the reader that the brightening considered here is not necessarily coincident with the GW chirp, and instead can be produced days to years after the merger when the cavity opened by the binary in the disc is refilled.

In conclusion, a significant fraction of our merger remnants (5 − 30%) can be detected by high-sensitivity instruments in X-rays and dominate the contaminating galactic emission, especially in the case of a merger-driven brightening. In the UV, galaxies strongly outshine the AGN, rendering almost the whole sample unobservable. We note that here we only consider photometry, but emission lines from the AGN could help disentangle the relative contributions if sufficiently luminous. Since in our model AGN can only be observed in the X-rays, henceforth we restrict our analysis to this band.

3.2.2. X-ray transients

Simply being above the detectability threshold is not a sufficient condition to identify a merger. Mergers also typically require a bright transient signature to ease their identification. In our model, we consider that the accretion rate can attain f_Edd = 1 (in a ‘brightening’) around the time of the merger, see Sect. 2.7.

We study the magnitude of merger-driven transient signals in our sample by comparing the pre-transient total flux with the brightening total flux for each source. Here, by total flux we mean the sum of fluxes from the remnant BH and its host galaxy. We recall that we assume that a merger-induced variation in the X-ray flux can be detected if either: (i) The source ‘appears’ – it is undetected before and detected at the brightening. The BH is observable at the brightening. (ii) The source is detected in both cases, but the flux changes significantly, by more than a factor of 2. This is probably an optimistic choice since AGN are intrinsically variable sources. Additionally, the BH has to be observable before or after. We denote mergers fulfilling any of these criteria as EM counterparts.

Transients are explored in Fig. 6. Pre-transient and brightening fluxes are connected by orange dashed or solid lines if there is an EM counterpart, corresponding observable transients of case (i) or (ii), respectively. We find that 4% of numerical mergers (37 out of 878) have an EM counterpart, according to any of our criteria. For delayed mergers, the fraction rises to almost 20%.

Fig. 6. X-ray transients in the numerical merger sample. Similar to the right panel of Fig. 5. Black dots and red crosses correspond to the X-ray post-merger BH fluxes for the pre-transient case (fiducial properties) and the brightening (which assumes fEdd = 1). We only show mergers for which either of the two configurations is observable. Pre-transient and brightening fluxes are connected by orange dashed lines if the pre-transient source (total flux from the BH and the galaxy) is undetected and the brightening source is detected and the BH is observable, and orange solid lines if both fluxes are detectable and the flux change is larger than a factor of 2. GW-undetectable remnant BHs are shown in purple (pre-transient) and orange (brightening).

As discussed above, a brightening can make observable a large number of AGN that would otherwise be too faint to be observed. A large fraction of EM counterparts (∼55%) are of type (i) in our notation, i.e. they ‘appear’ at the brightening. We note that many of these events are only marginally observable since our sample is dominated by low-flux sources. These mergers would in practice be hard to detect. Additionally, many of these mergers have very low (pre-transient) BH accretion rates or low mass ratios. The brightening luminosity could be much dimmer than estimated in many of these systems due to the low availability of gas. In this sense, our transient model is optimistic.

Around 45% of the sample is of type (ii), that is, the source is detectable both before and after the merger and the change in the flux at brightening is larger than 2. This sample is dominated by sources that are initially galaxy-dominated and then become AGN-dominated. In this case, the spectral shape of the source could also change at the brightening. In contrast, only in ∼10% of EM counterparts, the source is AGN-dominated before and at the brightening. Remnant BHs which are bright enough to be observable before the brightening, generally already accrete at f_Edd ∼ 1. Because of this, the flux difference generated by the brightening is low and such bright mergers are unlikely to have a detectable transient feature.

We note that it is likely that our model overestimates the pre-transient emission and thus underestimates the number of transients of type (ii) since the presence of the cavity can reduce the emission from the disc by more than a factor ∼2 (Bogdanović et al. 2022, and references therein). In our SED model, if the disc is truncated at 20GM/c² due to the cavity, the X-ray luminosity of all BHs with spin ≳0.5 (which encompasses almost all BHs observable at the brightening) decreases by a factor of > 2. Assuming a pre-merger decrease in the X-ray luminosity drops by a factor of 2 would lead instead to a factor of 2 increase in the number of mergers having a detectable transient. It is also possible the merger-induced luminosity burst exceeds the Eddington luminosity (Armitage & Natarajan 2002).

3.2.3. The population of X-ray-observable mergers

In the sections above, we identified several sub-samples of BH mergers that are more favourable for X-ray detection. These sub-samples do not reflect the properties of the global population of BH mergers. If future instruments used X-rays to detect BH mergers, they would be biased with respect to the global merger population. In this section, we quantify such biases by studying the differences between AGN-dominated mergers, observable mergers, EM counterparts, and the global population of BH mergers. We will use the numerical merger sample in order to account for the multi-messenger observability of BH mergers since the GW analysis was only performed for this sample.

Since more massive BHs are generally brighter, observable remnants are strongly biased towards massive mergers with M_• ≳ 10⁶ M_⊙ and M_* ≳ 10¹⁰ M_⊙, and sources with a high M_•/M_* ratio, which are more likely to be AGN-dominated. This results in AGN-dominated and observable mergers being on average over-massive with respect to the global BH merger population at fixed M_*, for 10^9.5 M_⊙ ≲ M_* ≲ 10^10.5 M_⊙ (Fig. 7). This effect is akin to the Lauer bias (Lauer et al. 2007) that causes AGN in a flux-limited sample to yield a relation between BHs and galaxy properties with BH masses above the ‘true’ relation for the full underlying population, especially at high redshift. The exception is observable and AGN-dominated remnants in M_*≳^10.5 M_⊙ galaxies, which consist of almost all bright high-mass mergers in our sample as they have systematically high accretion rates and M_•/M_* ratios and therefore are not biased.

Fig. 7. Mass of the remnant BH against the mass of the host galaxy for mergers occurring in the range z = 3.5 − 4. Blue crosses correspond to the overall population of numerical remnant BHs, crimson triangles to AGN-dominated mergers and green circles to observable mergers in the X-ray. Lines denote the geometric average of M• in bins of M*, while shaded regions denote the 1σ error in the average. We note that all observable mergers are AGN-dominated sources. AGN-dominated mergers are more massive, at fixed M* compared to all merger remnants: this is an effect of requiring the AGN to be bright, and brighter than the galaxy. This effect is stronger for observable mergers. A similar trend is found at all redshifts.

In Fig. 8, we compare the distributions in M_*, M_•, sSFR, f_Edd, redshift, and q of the global BH merger population, AGN-dominated sources, and observable remnant BHs.

Fig. 8. Distribution of M*, M•, sSFR, fEdd, redshift and q for the overall population of numerical BH mergers (in blue) and for the sub-sample of AGN-dominated mergers (in crimson), and observable mergers (in green) in the X-rays. The integrals of all distribution functions are normalised to unity. X-ray-observable merger remnants accrete faster and are more massive with respect to the general merger remnant population. They occur at lower redshifts and are hosted in more massive galaxies with lower sSFR.

We first look at the AGN-dominated sources in comparison to the global BH merger population. We find that, critically, very high accretion rates are needed for the remnant BH to be bright enough to dominate the galactic emission. Also, BHs hosted by galaxies with low XRB emission (i.e. sSFR) are more likely to be AGN-dominated. This biases the sample of AGN-dominated mergers towards lower sSFR galaxies at fixed M_* and towards high-mass galaxies, which generally have lower sSFR. Finally, high-mass BHs are only assembled at lower redshifts, and so AGN-dominated mergers tend to occur at lower redshifts.

Observable mergers have similar characteristics, further exacerbated by the flux requirement, which selects only BHs with M_• ≳ 10⁶ M_⊙. This leads to the selection of host galaxies with M_* ≳ 10¹⁰ M_⊙, low sSFR, and z ≲ 4.5. Observable remnants are likewise strongly biased towards highly-accreting BHs, with almost all remnants having f_Edd > 0.3. The mass ratios are small because numerical BH mergers involving a massive primary often have a much lighter secondary. This effect is much less pronounced for delayed mergers.

In summary, AGN-dominated and observable mergers are strongly biased towards highly-accreting BHs hosted in galaxies with low sSFR. Observable mergers are further biased towards high BH and host galaxy masses. We note that a significant fraction of the low-M_• AGN-dominated population would not be present if dynamical delays were taken into account since the delay times of BHs in low-mass galaxies are long, so the BHs would grow significantly during the delay or not coalesce before the end of the simulation.

In Fig. 9, we consider that the merger produces an Eddington-limited brightening around the time of the merger, increasing the luminosity. In this scenario, a large fraction of sources is AGN-dominated, and so AGN-dominated mergers trace the global merger population well. Observable BHs now include a larger number of lower mass mergers (M_• ∼ 10⁶ M_⊙ and M_* ∼ 10¹⁰ M_⊙) that would otherwise not be observable given their low fiducial accretion rates. An important caveat here is that the brightening luminosity might depend on the pre-transient accretion, which our model does not take into account.

Fig. 9. Similar to Fig. 8, but assuming a brightening luminosity such that fEdd = 1. In orange, we show the distribution of observable mergers according to any of the criteria defined in the text. The presence of a brightening can significantly increase the number of lower-mass X-ray-observable mergers. EM counterparts are biased in a similar way to X-ray-observable remnants at the brightening.

The population of EM counterparts is similar to that of observable brightenings. Nevertheless, some of the brightest, most massive BHs are not detected as EM counterparts. This is because high-mass BHs accrete at high rates before the brightening, so the assumed f_Edd = 1 brightening does not increase the flux significantly (see Sect. 3.2.2). We note that, as discussed in Sect. 3.2.2, this population of bright massive mergers could be observable if the presence of the inner disc cavity significantly decreases the pre-transient luminosity.

The population of X-ray-observable mergers can have on average worse GW detectability and parameter estimation than the global merger sample. X-ray-observable mergers tend to be biased towards high-mass mergers, which tend to have low mass ratios. Such high-mass, low-q mergers tend to be harder to detect with LISA. In Fig. 10, we study the GW sky localisation error of the electromagnetically detectable mergers in our numerical merger sample. The GW sky localisation of EM counterparts is similar to that of the global numerical merger sample. 31% of X-ray EM counterparts are localised with 90% confidence accuracy smaller than 10 deg², which is comparable with the value of 37% for all mergers. This fraction is lower (15%, 6 out of 41) for X-ray-observable remnants tends since they are more strongly biased towards high-mass low-q mergers which are poorly detected. We find that most X-ray-bright mergers can in general also be detected by LISA, with a S/N distribution comparable to that of the global merger population. Of the 46 X-ray-observable mergers, 41 are also GW-detected.

Fig. 10. LISA 90%-confidence GW sky localisation error of EM-observable mergers. The top panel shows the distribution of all numerical mergers (in blue), X-ray-observable mergers (in green), and radio-observable mergers with SKA in the pessimistic model (in pink). The bottom panel shows the distribution of X-ray EM counterparts (in orange) and radio EM counterparts (in yellow).

In conclusion, the population of AGN-dominated sources, observable sources, and EM counterparts are biased tracers of the underlying merger population. Observable BHs and mergers are more massive, inhabit more massive and less star-forming galaxies, accrete at higher rates, and occur at lower redshifts. They also tend to be overmassive with respect to the M_•-M_* relation.

Finally, we remind the reader that our discussion is limited to high-z events since the simulation stops at z ∼ 3.5. At lower redshift, a larger fraction of the merger population could be detectable since the sources will be closer to the observer. Furthermore, we have not included delayed mergers here, whose lower redshift and mass ratios closer to unity would improve the S/N and sky localisation.

3.3. Radio

3.3.1. Radio detectability

The radio observability of merger remnants jets at ν = 2 GHz (observer-frame) is explored in Fig. 11. We recall that we consider two models: a lower limit model for the core radio emission based on the fundamental plane (following Gültekin et al. 2009), and an upper limit model for the total radio luminosity based on the theoretical model in Meier (2001). We also recall that we do not explicitly consider the contamination due to the galactic radio emission in our analysis since we cannot quantify this reliably (see Sect. 2.8). Some of the sources which exceed the instrumental sensitivity threshold could be outshone by their hosts if M_• ≲ 10^7.5 M_⊙. In this section, we disregard this effect and denote remnant BHs as observable if their flux exceeds the instrumental sensitivity and for numerical mergers if they are also detected by LISA.

Fig. 11. Observer-frame spectral flux density at ν = 2 GHz of numerical merger remnants against the remnant BH mass, for numerical mergers (left panel) and delayed mergers (right panel). For each sample, the ‘lower limit’ model (based on the fundamental plane, Gültekin et al. 2009) is shown with dark red circles, and the ‘upper limit’ model (based on the theoretical model in Meier 2001) is shown with dark green crosses. The dashed horizontal lines represent the 5σ sensitivity thresholds for the ngVLA (Carilli et al. 2015), in black, for the SKA1-MID, in grey, and for the SKA (Prandoni & Seymour 2015), in purple. These thresholds are calculated assuming 9 h exposure times and that the observed sources are not resolved. The legend indicates the fraction of observable events, which lie above the sensitivity threshold of each instrument and are observable with LISA, if applicable. Mergers that are undetected with LISA are shown with bigger and lighter markers.

We consider the detectability by future surveys with ngVLA for which we consider a sensitivity threshold 1.5 μJy (Carilli et al. 2015) and SKA1-MID with a threshold of 4.2 μJy (Prandoni & Seymour 2015). The thresholds are calculated assuming 9 h exposure observations at 2 GHz and that sources are not resolved. We also consider that the full SKA could go deeper, and assume a sensitivity of 0.4 μJy. These detectability limits are denoted by grey, black, and purple dashed lines in Fig. 11 respectively.

We find that a fraction of mergers can be detected in the radio with future instruments, although this fraction depends strongly on the model assumed for the radio luminosity and on the instrument’s sensitivity. For the pessimistic empirical model for core emission, only the most massive BHs with masses M_• ≳ 10^7.5 M_⊙ can be detected. The fraction of observable mergers is in the range 1 − 3%, which is lower than the fraction we found for X-rays in Sect. 3.2.1. For the optimistic theoretical model for the radio emission, BHs with mass M_• ≳ 10^6.5 M_⊙ are generally above the flux limit. The fraction of observable mergers rises significantly to 4 − 30%, consistent with the fractions found for the X-rays detectability. The different predictions for our two models stem from the fact that the pessimistic model estimates only the core emission while our optimistic model estimates the total emission, although it is also possible that the pessimistic model underestimates the radio luminosity for highly accreting AGN (Gültekin et al. 2022).

3.3.2. Radio transients

In radio, we consider two possible types of transients. Firstly, analogous to our model for X-ray transients, we consider a merger brightening caused by a sharp increase in the accretion rate. Secondly, we consider a flare, as found in some simulations (see Sect. 2.7 for more details). As for the X-ray transients, we consider three possibilities for detecting an EM counterpart, (i) The source is detected only after the merger, (ii) The source is detected before and after, but the flux changes by more than a factor of 2. As in the previous section, we do not take into consideration the galactic contribution to the flux. Henceforth, we take the fundamental plane pessimistic model as the fiducial radio model for this analysis since it models only the core luminosity, which is expected to show stronger flux variations on timescales of weeks to months, and we optimistically take SKA as the fiducial instrument. Since the GW analysis was only performed on numerical mergers, we use numerical mergers as the reference sample.

In Fig. 12 we show results for the brightening (top panel) and the flare (bottom panel). We find that observing a radio EM counterpart is unlikely – less than 1% of mergers have a detectable merger-induced transient. ∼80% of EM counterparts are observable due to a post-merger brightening, i.e. an increase in the accretion rate to f_Edd = 1. Most potentially observable brightenings already accrete at very high rates before the merger and so the merger-induced change in luminosity is small, which means such mergers will not be recognisable as radio transients. Transients EM counterparts are rarer in the radio than in the X-ray as at fixed BH mass L_R, aft/L_R = (L_X, aft/L_X)^0.67 (Eq. (14)), where the suffix ‘aft’ denotes the brightening luminosity, the flux variation due to the brightening in radio is smaller than in X-rays. Only one EM counterpart (∼20% of all EM counterparts, bottom panel of Fig. 12) is observable due to a flare. The amplitude of the flare is strongly dependent on the mass ratio, and thus only very bright and massive mergers with nearly equal mass ratios, which are very uncommon in our sample, can have an observable flare. In general, the number of EM counterparts is small also because the number of observable BHs is small.

Fig. 12. Observability of a transient generated by a fEdd = 1 brightening (top panel) or a merger-induced flare (bottom panel) in the numerical merger sample. Black dots and red crosses correspond to the post-merger BH radio fluxes for the fiducial and the rebrightening scenario, respectively. The pre-merger primary BH radio flux is shown in blue triangles, and the flare flux is shown in green squares. We only show mergers for which either the fiducial or the merger-enhanced (brightening or flare) flux is observable. Fiducial and merger-enhanced fluxes are connected by orange dashed lines if the fiducial flux is undetected and the enhanced flux is detected, or orange solid lines if the flux is detected in both cases and the flux change is larger than a factor of 2. For simplicity, we only consider the detectability of the pessimistic model by the SKA (assuming a sensitivity of 0.4 μJy, purple horizontal line), but the ngVLA sensitivity is also shown for reference (1.5 μJy, black horizontal line). GW-undetected events are shown with bigger, brighter markers. With our modelling, most of the merger-induced flux changes are too small to be detectable in radio.

Finally, we note that these estimates using the pessimistic model based on the fundamental plane are likely lower limits since the fundamental plane can underestimate the radio luminosity of highly accreting BHs (Gültekin et al. 2022).

3.3.3. The population of radio-observable mergers

The population of radio-observable remnant BHs is also a biased tracer of the underlying remnant BH population. In order to quantify this, in Fig. 13 we compare the properties of observable mergers and mergers with transient EM counterparts to the global BH merger population. We also show the population of radio-observable mergers which are also X-ray-observable. Again we use the pessimistic model and SKA as the fiducial model and instrument. For a total of 878 numerical mergers, there are 21 radio-observable mergers and 5 EM counterparts. The biases of the radio-observable population are qualitatively identical to those of the X-ray observable population – observable BHs have on average significantly higher M_•, M_* and f_Edd and lower sSFR and redshift. However, given that a smaller amount of mergers are observable for our pessimistic model, the sample is more strongly biased towards high BH and galaxy masses (M_• ≳ 10⁷ M_⊙, M_* ≳ 10¹⁰ M_⊙) and low redshifts (z ≲ 4). Radio-observable mergers are also overmassive with respect to their galaxies. We find that all radio-observable mergers in our sample are also X-ray-observable since X-ray bright BHs tend to also be bright in the radio.

Fig. 13. Distribution of M*, M•, sSFR, fEdd, redshift, and q for the overall population of numerical BH mergers (in blue), for radio-observable mergers (AGN-dominated mergers where FX, • > FX, lim, in green) and for EM counterparts (in orange). All histograms are normalised to unity. Observable BH mergers are hosted in more massive galaxies, have lower sSFR and higher fEdd with respect to the general BH-merging population. All radio-observable mergers are also X-ray-observable.

As discussed in the previous section, the number of EM counterparts is very small since neither the brightening nor the flare is able to sufficiently change the flux for the brightest BHs and produce an observable transient. The population of EM counterparts is similar to that of observable mergers, although the number of these mergers is too low to carry out a meaningful statistical analysis.

Since radio-observable mergers are strongly biased towards high-mass low-q mergers in our pessimistic model, they have poorer GW parameter estimation with LISA. As shown in Fig 10, few radio-observable remnants or radio EM counterparts have a 90%-confidence error lower than 10 deg². Despite this, radio-observable remnants and radio EM counterparts tend to be detectable with LISA, although GW-undetected mergers constitute a fraction of ∼20% (5 out of 26 radio-observable mergers). Finally, we refer the reader again to the note at the end of Sect. 3.2.3, which applies also to radio sources. As before, we note that the results presented in this section are likely lower limits to the radio observability of merger remnants since our pessimistic model is likewise a lower limit on the luminosity of highly accreting BHs (Gültekin et al. 2022).

4. Comparison with previous work

Previous work by Tamanini et al. (2016) and Mangiagli et al. (2022) also studied the possibility of joint (EM and GW) multi-messenger detections of BH mergers. They use a BH population synthesised from a semi-analytic model of galaxy formation. The X-ray AGN fluxes estimated by Mangiagli et al. (2022) are similar to those presented in this work. The fraction of detected radio sources found in the present work is significantly lower than the values reported in Tamanini et al. (2016) and Mangiagli et al. (2022), despite the optimistic model presented here being conceptually similar. Our luminosities can be more than 4 orders of magnitude smaller, for two reasons: Firstly, for our optimistic model we assume a conversion factor from jet power to radio luminosity of 10⁻², instead of assuming full or very efficient conversion, and take into account the fact that the radiation is distributed across a wide synchrotron spectrum. Secondly, we use a different model for the flare emission, which predicts smaller luminosities. Lops et al. (2023) used similar methods as Tamanini et al. (2016) and Mangiagli et al. (2022) for the synthesis but focus on the GW localisation of BH mergers, by studying the galaxy fields in LISA error-boxes.

Other previous studies (e.g., Kelley et al. 2019; Krolik et al. 2019) have rather focused on BH binary signatures, such as a periodic modulation of the BH luminosity or spectral features. Many of these signatures are expected already well before the BH coalescence and in only few cases they can be observed concurrently with a GW detection with LISA, although they can provide valuable information to break degeneracies and compare the speed of photons and gravitons (Haiman 2017).

The main novelty of the present study with respect to previous work is that we use a hydrodynamical simulation that follows self-consistently the evolution of BHs and their environments and obtain a realistic BH merger population. This also allows us to calculate consistently environmental parameters such as gas and dust obscuration. We also introduce physically motivated models for the AGN SED, radio emission, and subgrid obscuration. Further, we study several important observational effects: (i) the contamination from the host galaxy emission (ii) the observability of a possible merger transient signal occurring around the time of the merger (iii) the observational biases of the EM-observable population.

5. Conclusions

In this work, we have presented a comprehensive study of BH mergers in OBELISK, a cosmological hydrodynamical simulation following the evolution of a protocluster down to redshift z ∼ 3.5. Building on Dong-Páez et al. (2023), which studied the properties of the population of BH mergers and compared it with the underlying global population of main BHs, we have performed a multi-messenger analysis of the detectability of BH mergers, in order to forecast their detectability and assess the possible observational biases. We performed a GW and EM analysis of numerical BH mergers at the resolution of the simulation and an EM analysis of delayed BH mergers which consider post-processed dynamical delays below the simulation resolution. We summarise our results below:

• Most of the numerical merger sample (∼99%) can be detected by LISA, generally with very high S/N (Fig. 4). Only a small sample of high-mass low-q mergers and low-mass high-z mergers are undetected. The intrinsic binary parameters, such as the BH masses, spins, and redshift, can generally be measured with high accuracy. Only ∼37% of these high-z BH mergers can be localised in the sky with 90%-confidence error better than 10 deg².

• In UV, remnant BHs are significantly fainter than their host galaxies, which are actively star-forming at the redshifts considered (z > 3.5).

• In X-rays, 5 − 15% of remnant BHs are bright enough to be detectable by future instruments while dominating over their host galaxy’s emission (Fig. 5). If a merger-induced brightening increases the BH accretion rate to the Eddington rate, up to 10 − 30% could become observable and in some cases be identifiable as a transient EM counterpart (Fig. 6).

• X-ray observable BH remnants tend to accrete near the Eddington limit, have higher M_* (≳10¹⁰ M_⊙), and M_• (≳10⁶ M_⊙), occur at lower redshift, in lower sSFR galaxies (Figs. 8 and 9), and are overmassive at fixed galaxy mass with respect to the global population for M_* ≲ 10^10.5 M_⊙ (Fig. 7)

• In radio between ∼1% and 30% of mergers can be detected by future instruments (Fig. 11). A merger-induced increase of the core radio luminosity is found to lead to only a small number of transients detectable as EM counterparts (Fig. 12).

• The population of radio-observable BH merger remnants differ with respect to the full merger population. The biases with respect to the full merger population are qualitatively analogous to the X-ray observable sample (Fig. 13). Most radio-observable mergers are also X-ray observable.

Overall, we found that the number of EM counterparts is currently limited by LISA’s ability to localise the systems in the sky. It is worth noting that an additional LISA-like detector would dramatically change the situation, leading to a sky localisation improvement of two orders of magnitude (Ruan et al. 2020b) in which case we would instead be limited by the sensitivity of EM telescopes and by the presence or absence of a tell-tale transient EM sign of BH mergers.

The expectation is that many if not most LISA BH mergers will have z ≳ 2, therefore our work considers some of the most plausible types of sources, but high-redshift sources are by definition fainter than low-redshift sources at fixed luminosity. At lower redshifts, prospects for looking for EM counterparts are brighter (Lops et al. 2023). In future work, we will explore the detectability of BH mergers in OBELISK at additional wavelengths.

Acknowledgments

We thank Geoffrey Bicknell, Sera Markoff, Stanislav Babak, and Alexander Wagner for stimulating discussions, and Sylvain Marsat for allowing usage of the lisabeta code for this paper. We thank the anonymous referee whose comments and suggestions helped improve the manuscript. MV, YD, NW, and SV acknowledge funding from the French National Research Agency (grant ANR-21-CE31-0026, project MBH_waves). MT acknowledges support from the NWO grant 0.16.VIDI.189.162 (‘ODIN’). AM acknowledges support from the postdoctoral fellowships of IN2P3 (CNRS). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101066346 (MASSIVEBAYES). This work has received funding from the Centre National d’Etudes Spatiales. This work has made use of the Horizon Cluster hosted by Institut d’Astrophysique de Paris; we thank Stéphane Rouberol for running smoothly this cluster for us. We acknowledge PRACE for awarding us access to Joliot Curie at GENCI@CEA, France, which was used to run most of the simulations presented in this work. Numerical computations were partly performed on the DANTE platform, APC, France. Additionally, this work was granted access to the HPC resources of CINES under allocations A0040406955 and A0040407637 made by GENCI.

References

All Figures

	Fig. 1. Gas density projection of a region of the OBELISK simulation at redshift 4.02. BHs with M• > 106 M⊙ are denoted by black star symbols.
In the text

Fig. 2. Spectral energy distribution of a merger remnant with mass 1.3 × 108 M⊙, fEdd = 0.51 and spin 0.90, residing in a galaxy with mass 1.4 × 1011 M⊙ at redshift 3.55. The top panel depicts the (unobscured) rest-frame luminosities while the bottom panel depicts the (obscured) observer-frame flux. The SED of the remnant BH is shown in the solid black and red lines for the fiducial and the brightening scenario, while the approximate UV and X-ray galaxy SED is shown in the blue dashed lines. The green arrows indicate approximately the instrumental sensitivities assumed in radio at 2 GHz (assuming a sensitivity of 0.4 μjy), in UV at 1500 Å (mUV, lim = 31.3), and in the 0.5 − 2 keV band (FX, lim = 10−17 erg s−1).

In the text

Fig. 3. Distribution of the optical depth against the host galaxy stellar mass for remnant BHs in our numerical merger sample. Top panel: rest-frame UV; middle panel: observer-frame 0.5 − 2 keV; bottom panel: observer-frame 2 − 10 keV. The colour denotes the number of BHs in each bin, on a logarithmic scale. Optical depths are very high in the UV, while in X-rays the optical depths are not extreme.

In the text

Fig. 4. Parameter distribution for the mock LISA data analysis. The top-left panel shows the signal-to-noise ratio (S/N) distributions. Gravitational wave events are assumed to be detected if S/N > 10, indicated by the vertical dashed line. The 90%-confidence error in the sky localisation is shown in the bottom-left panel. The dashed line delimits the shaded region, which corresponds to the events with sky localisation poorer than 10 deg2. The relative 90%-confidence error distributions for z, Mchirp and q, and for χ1 and χ2, are shown in the right panels, top and bottom respectively. Most parameters are recovered with small uncertainties. The majority of mergers have sky localisation worse than 10 deg2 at the time of the merger.

In the text

Fig. 5. Black-hole vs galaxy flux, both corrected for absorption. Left column: rest-frame UV magnitude at ν = 1500 Å; right column: X-ray 0.5 − 10 keV band. Top row: numerical mergers; bottom row: delayed mergers. Black dots: accretion rate measured in the simulation; red crosses: merger-induced brightening with fEdd = 1. The grey region corresponds to galaxies dominating the emission over their BHs. The black dashed line indicates the boundary where the luminosity and galaxy and BH are equal. The hatched region corresponds to galaxies and BHs below the minimum flux needed for detection, which we set to mUV, lim = 31.3 (UV) and FX, lim = 10−17 erg cm−2 s−1 (X-rays). In the top left corner, we indicate the fraction of BHs with higher flux than their host galaxies and that of BHs above the observational flux limit. We note that a significant fraction of sources are very faint and lie beyond the limits of the plot, especially in the UV. Mergers that are GW undetected by LISA are shown with bigger and lighter markers, in purple and orange, for the fiducial and brightening cases. In UV all merging BHs are outshone by their host galaxy, while in X-rays a fraction of merging BHs are both brighter than the host galaxy and detectable by sensitive X-ray telescopes.

In the text

Fig. 6. X-ray transients in the numerical merger sample. Similar to the right panel of Fig. 5. Black dots and red crosses correspond to the X-ray post-merger BH fluxes for the pre-transient case (fiducial properties) and the brightening (which assumes fEdd = 1). We only show mergers for which either of the two configurations is observable. Pre-transient and brightening fluxes are connected by orange dashed lines if the pre-transient source (total flux from the BH and the galaxy) is undetected and the brightening source is detected and the BH is observable, and orange solid lines if both fluxes are detectable and the flux change is larger than a factor of 2. GW-undetectable remnant BHs are shown in purple (pre-transient) and orange (brightening).

In the text

Fig. 7. Mass of the remnant BH against the mass of the host galaxy for mergers occurring in the range z = 3.5 − 4. Blue crosses correspond to the overall population of numerical remnant BHs, crimson triangles to AGN-dominated mergers and green circles to observable mergers in the X-ray. Lines denote the geometric average of M• in bins of M*, while shaded regions denote the 1σ error in the average. We note that all observable mergers are AGN-dominated sources. AGN-dominated mergers are more massive, at fixed M* compared to all merger remnants: this is an effect of requiring the AGN to be bright, and brighter than the galaxy. This effect is stronger for observable mergers. A similar trend is found at all redshifts.

In the text

Fig. 8. Distribution of M*, M•, sSFR, fEdd, redshift and q for the overall population of numerical BH mergers (in blue) and for the sub-sample of AGN-dominated mergers (in crimson), and observable mergers (in green) in the X-rays. The integrals of all distribution functions are normalised to unity. X-ray-observable merger remnants accrete faster and are more massive with respect to the general merger remnant population. They occur at lower redshifts and are hosted in more massive galaxies with lower sSFR.

In the text

	Fig. 9. Similar to Fig. 8, but assuming a brightening luminosity such that fEdd = 1. In orange, we show the distribution of observable mergers according to any of the criteria defined in the text. The presence of a brightening can significantly increase the number of lower-mass X-ray-observable mergers. EM counterparts are biased in a similar way to X-ray-observable remnants at the brightening.
In the text

	Fig. 10. LISA 90%-confidence GW sky localisation error of EM-observable mergers. The top panel shows the distribution of all numerical mergers (in blue), X-ray-observable mergers (in green), and radio-observable mergers with SKA in the pessimistic model (in pink). The bottom panel shows the distribution of X-ray EM counterparts (in orange) and radio EM counterparts (in yellow).
In the text

Fig. 11. Observer-frame spectral flux density at ν = 2 GHz of numerical merger remnants against the remnant BH mass, for numerical mergers (left panel) and delayed mergers (right panel). For each sample, the ‘lower limit’ model (based on the fundamental plane, Gültekin et al. 2009) is shown with dark red circles, and the ‘upper limit’ model (based on the theoretical model in Meier 2001) is shown with dark green crosses. The dashed horizontal lines represent the 5σ sensitivity thresholds for the ngVLA (Carilli et al. 2015), in black, for the SKA1-MID, in grey, and for the SKA (Prandoni & Seymour 2015), in purple. These thresholds are calculated assuming 9 h exposure times and that the observed sources are not resolved. The legend indicates the fraction of observable events, which lie above the sensitivity threshold of each instrument and are observable with LISA, if applicable. Mergers that are undetected with LISA are shown with bigger and lighter markers.

In the text

Fig. 12. Observability of a transient generated by a fEdd = 1 brightening (top panel) or a merger-induced flare (bottom panel) in the numerical merger sample. Black dots and red crosses correspond to the post-merger BH radio fluxes for the fiducial and the rebrightening scenario, respectively. The pre-merger primary BH radio flux is shown in blue triangles, and the flare flux is shown in green squares. We only show mergers for which either the fiducial or the merger-enhanced (brightening or flare) flux is observable. Fiducial and merger-enhanced fluxes are connected by orange dashed lines if the fiducial flux is undetected and the enhanced flux is detected, or orange solid lines if the flux is detected in both cases and the flux change is larger than a factor of 2. For simplicity, we only consider the detectability of the pessimistic model by the SKA (assuming a sensitivity of 0.4 μJy, purple horizontal line), but the ngVLA sensitivity is also shown for reference (1.5 μJy, black horizontal line). GW-undetected events are shown with bigger, brighter markers. With our modelling, most of the merger-induced flux changes are too small to be detectable in radio.

In the text

Fig. 13. Distribution of M*, M•, sSFR, fEdd, redshift, and q for the overall population of numerical BH mergers (in blue), for radio-observable mergers (AGN-dominated mergers where FX, • > FX, lim, in green) and for EM counterparts (in orange). All histograms are normalised to unity. Observable BH mergers are hosted in more massive galaxies, have lower sSFR and higher fEdd with respect to the general BH-merging population. All radio-observable mergers are also X-ray-observable.

In the text