Detectability of a spatial correlation between stellar-mass black hole mergers and Active Galactic Nuclei in the Local Universe

The origin of the Binary Black Hole (BBH) mergers detected through Gravitational Waves (GWs) by the LIGO-Virgo-KAGRA (LVK) collaboration remains debated. One fundamental reason is our ignorance of their host environment, as the typical size of an event's localization volume can easily contain thousands of galaxies. A strategy around this is to exploit statistical approaches to assess the spatial correlation between these mergers and astrophysically motivated host galaxy types, such as Active Galactic Nuclei (AGN). We use a Likelihood ratio method to infer the degree of GW-AGN connection out to $z=0.2$. We simulate BBH mergers whose components' masses are sampled from a realistic distribution of the underlying population of Black Holes (BHs). Localization volumes for these events are calculated assuming two different interferometric network configurations. These correspond to the configuration of the third (O3) and of the upcoming fourth (O4) LVK observing runs. We conclude that the 13 BBH mergers detected during the third observing run at $z\leq0.2$ are not enough to reject with a \(3\sigma\) significance the hypothesis according to which there is no connection between GW and AGN more luminous than $\approx 10^{44.3}\rm{erg}\ \rm{s}^{-1}$, that have number density higher than \(10^{-4.75}\textrm{Mpc}^{-3}\). However, 13 detections are enough to reject this no-connection hypothesis when rarer categories of AGN are considered, with bolometric luminosities greater than $\approx 10^{45.5}\rm{erg}\ \rm{s}^{-1}$. We estimate that O4 results will potentially allow us to test fractional contributions to the total BBH merger population from AGN of any luminosity higher than \(80\%\).


INTRODUCTION
Since the first direct GW detection has been announced (Abbott et al. 2016), the two interferometers of Advanced LIGO (LIGO Scientific Collaboration et al. 2015) and the one of Advanced Virgo (Acernese et al. 2015) have measured the signal coming from tens of compact objects mergers in three observing runs (Abbott et al. 2019(Abbott et al. , 2021aThe LIGO Scientific Collaboration et al. 2021a,b). Thanks to improved sensitivities and the addition of a fourth detector, KA-GRA (Somiya 2012;Aso et al. 2013), this number will grow in the upcoming years (Abbott et al. 2018).
Different formation pathways for these merging BBHs have been proposed (Mapelli 2021). They might arise from the evolution of isolated stellar binary systems (Dominik et al. 2012;Belczynski et al. 2016;Spera et al. 2019), or form in dense environments, in which dynamical interactions can efficiently drive binaries of compact objects towards the merger (Stone 2017;Rodriguez & Loeb 2018;Antonini et al. 2019;Gerosa & Berti 2019;Rodriguez et al. 2021;Rizzuto et al. 2021). One particular example of such environments can be the accretion disk around Super Massive Black Holes (SMBHs) in Active Galactic Nuclei (Bartos et al. 2017b;Stone et al. 2017;McKernan et al. 2018;Ford et al. 2019;Samsing et al. 2020;Gayathri et al. 2021). It has been shown that in such an environment, compact ob-★ E-mail: veronesi@strw.leidenuniv.nl jects can migrate towards a radius close to the one of the innermost stable circular orbit, and there be trapped for the remaining AGN lifetime (Peng & Chen 2021). The large number density of compact objects and the high escape velocity in that region can facilitate the occurrence of hierarchical mergers (i.e. mergers in which at least one of the two components is the remnant of a previous merger) (Yang et al. 2019;Gerosa & Fishbach 2021;Wang et al. 2021b). The mass of the remnants of hierarchical mergers can be higher than the lower bound of the Pair Instability mass gap predicted by stellar formation models (Farmer et al. 2019;Woosley & Heger 2021). This formation pathway has therefore the theoretical advantage of being able to explain the non-vanishing merger rate inferred for binaries with components heavier than 50M (The LIGO Scientific Collaboration et al. 2021c).
There are potentially several ways to address the formation pathways' open question and in particular to assess a plausible connection between GW events (BBH mergers in particular) and AGN. The most straightforward would be to directly detect an ElectroMagnetic (EM) counterpart in coincidence with the GW event. This might be possible in dense environments like the accretions disks of AGN (McKernan et al. 2019;Wang et al. 2021a), and such a counterpart might have already been observed (Graham et al. 2020) (However, see also Ashton et al. 2021). The typical localization volumes of GW events make their association with an EM counterpart challenging. The interferometers currently operating are in fact only able to as-sociate to GW detections comoving volumes that can easily contain thousands of different galaxies. Similarly to what happens in the case of the emission of an EM counterpart, a companion GW signal can be originated from the same source of a detected event in the case of mergers happening near a SMBH. These events could therefore be identified by the independent detection of an associated gravitational echo (Kocsis 2013;Gondán & Kocsis 2021).
Another way to infer the origin of the detected events is by statistically comparing the measured source population properties, such as mass and spin distributions, with model expectations. This kind of analysis has been done for several potential host environments, including AGN (McKernan et al. 2020;Gayathri et al. 2021;Tagawa et al. 2021;Wang et al. 2021b;Li 2022). While the expected distribution of spin parameters is still debated, all these analyses conclude that heavy (≥ 50 ) stellar-mass BHs are expected to be generated through the AGN formation channel.
Finally, the increasing number of detections allows us to exploit statistical approaches to explore the spatial correlation between GW events and specific types of host environments. These approaches can overcome the big challenge of large localization volumes. Bartos et al. (2017a) proposed a statistical likelihood-ratio-based method to find out how many GW detections would be needed to establish which fraction of BBH mergers detected through GWs happened in an AGN. This earlier work was based on the GW localization volume distribution expected for detections performed by the LIGO-Virgo network at design sensitivity and assuming only mergers of pairs of 10M BHs. In this work, we present an analysis based on the same method, although we use simulated GW detections constructed from the latest results on the inference of the underlying BBH component masses' distribution. To simulate these detections we employ detectors' sensitivities representative of the third observing run of the LIGO and Virgo interferometers, as well as those expected to characterize the fourth one, when KAGRA will join the network. This paper is organized as follows: In Section 2 we provide an overview of all the steps of the analysis, with details in the following subsections. How we simulated the GW detections the localization volumes of which are needed in the statistical analysis is described in subsection 2.1, while in subsection 2.2 we present how this statistical investigation works. The results of our works are presented ins section 3. Finally, in Section 4 we draw conclusions and discuss the next steps to improve our it to observed data.

METHOD
To investigate the spatial correlation between Gravitational Waves 90% Credibility Level localization volumes (thereafter "V90") and the positions of AGN in the local Universe, we first build two catalogues of simulated GW detections anchored in current observations. For the first, we simulate the response of the detector network active during O3. For the second catalogue, we use the same synthetic population of BBHs, and we simulate their detection by the interferometric network configuration expected for O4, which includes also KAGRA. To create the simulated detections we first sample the joint probability distribution of the binary mass ratio = 2 / 1 and primary component's mass 1 ; which is, by definition, greater than the mass of the secondary one, 2 . We then sample the spin distribution for each binary component, the distribution for the inclination of the orbital plane with respect to the line of sight, and for the luminosity distance between the position of the event and the detectors. The assumed distributions, as well as the configurations and the detector sensitivity curves used in our simulations, are described in section 2.1.
Once the mock observations have been simulated, we evaluate V90 for all the detections using BAYESTAR (Singer & Price 2016), a sky localization algorithm able to perform in a few seconds a Bayesian, non-Markov Chain Monte Carlo analysis.
We then use the newly created distribution for V90 to sample a set of comoving volumes that are then exploited in an algorithm based on the likelihood-ratio method described in Bartos et al. (2017a). This algorithm crossmatches the positions of the GW localization volumes with the positions of AGN in the local Universe, which are assumed to be isotropically distributed in comoving volume. The final output of this algorithm (described in detail in 2.2) is the number of GW detections needed to test the hypothesis according to which a certain fraction ( agn ) of the detected BBH mergers happened in an AGN; having the chance of rejecting the no-correlation hypothesis (none of the detected BBH mergers happened in an AGN) with a given confidence.

Simulation of GW detections
A distribution of V90 is required by our statistical method. We obtain such distribution by simulating several realistic GW detections for both O3 and O4 configurations. We describe the details of these simulations in the following sections.

Source population
Our simulated GW events are derived from the population analysis based on the latest results of the LVK Collaboration. We assume for 1 the P L + P analytical model presented in Abbott et al. (2021b) and we simultaneously sample values of 1 and from their joint posterior probability distribution. This distribution has been obtained through the standard hierarchical bayesian analysis presented in The LIGO Scientific Collaboration et al. (2021c) and posterior samples are publicly available. The secondary component's mass is then calculated as 2 = 1 . We use the same mass distribution to simulate BBH mergers irrespective of them happening in an AGN. This is done to maintain our estimate conservative and model-independent. We, therefore, neglect the effects of the hypotheses according to which GW events originated in dense environments are more likely to involve higher-mass BHs with respect to the ones that originated from an isolated binary system.
For simplicity, we assume for all the BHs the spin direction to be aligned with respect to the binary orbital angular momentum, and a uniform spin magnitude distribution between 0 and 1. The distribution of V90 is not expected to be significantly affected by such an assumption.
The simulated binaries are uniformly distributed in comoving volume and their inclinations are sampled according to a uniform distribution over arccos( ). The cosmological parameters we assume during our analysis are the ones inferred from the Planck Cosmic Microwave Background observations (Planck Collaboration et al. 2020).

The network of detectors
Next, we simulate the network of detectors. The whole analysis presented in this work is done for two different settings: the first one aims at reproducing the V90 distribution for O3, while the second one aims to forecast the distribution of the detected volumes expected during O4. In both cases, we assume a duty cycle of 78% for all the different detectors individually, and we keep a network Signal-to-Noise Ratio (SNR) threshold of 8, adding a Gaussian measurement error to the SNR and requiring that at least two detectors contribute to the network SNR with an individual SNR ≥ 4. The signals of the injected events are then compared with the detectors' noise in the 10 − 5000Hz frequency range. We use an IMRPhenomD waveform type Khan et al. 2016) to model the injections.
To reproduce the volume distribution of the events measured during O3, we model a network of three detectors: LIGO Hanford, LIGO Livingstone, and Virgo, using the sensitivities characterized by the following IDs: LIGOM L S P1200087 for the two LIGO interferometers, and A VM L S P1200087 for Virgo interferometer.
For the O4 predictions, we add a fourth KAGRA-like interferometer, and we change the sensitivity curves of each detector. Specifically, we use LIGOA VO4T1800545 for LIGO and Virgo detectors, and LIGOKAGRA80M T1800545 for KAGRA.

Evaluation of 90% CL localization volumes
For O3 (O4), out of the 200k (100k) injections, 663 (1737) have a SNR higher than the threshold. Out of these simulated mergers whose signals exceed the SNR threshold (hereafter referred to as detections), 274 for O3 and 317 for O4 have a measured value for the luminosity distance that corresponds to ≤ 0.2. We evaluate the value of V90 for each of these low-redshift events using the BAYESTAR algorithm (Singer & Price 2016). For these close events, we show the cumulative distribution of V90 in Figure 1. The blue and green histograms are detections simulated for O3 and O4, respectively. The top axis shows the expectation value of the number of AGN within the corresponding localization volume, assuming a uniform number density of AGN equal to agn = 10 −4.75 Mpc −3 . The same value for this parameter was used in Bartos et al. (2017a) and Corley et al. (2019). This number density corresponds to AGN with a bolometric luminosity higher than ≈ 10 44.3 erg s −1 in the local Universe. This value for the minimum bolometric luminosity for AGN at a specific number density has been obtained by integrating the double power law that represents the AGN L F in (Hopkins et al. 2007), using the best fit values for = 0.1. This holds for all the values of bolometric luminosities mentioned hereafter.
As a sanity check, we verify that our sample of V90 from O3 simulations and the values of V90 for the 13 observations of O3 with redshift ≤ 0.2 are compatible with a single common distribution. We do this with a 2 samples Kolmogorov-Smirnov test. We find that the hypothesis according to which the two samples come from the same distribution cannot be rejected (p-value ≈ 0.39).

Minimum number of GW detections to test the AGN origin
We consider a Universe where a fraction of GW events agn originate in an AGN-type of galaxy. Our goal is to calculate how many GW detections are needed to infer this AGN-BBH mergers connection; more precisely, the minimum number 3 gw of GW detections below = 0.2 needed to reject with a 3 significance the hypothesis of noconnection between detected BBH mergers and AGN. We evaluate 3 gw as a function of the fraction agn of GW events originated from an AGN. We calculate such a number by investigating the spatial correlation between AGN positions (assumed to be uniformly distributed in comoving volume) and the localization volumes of simulated GW detections, starting from the statistical approach presented in Bartos et al. (2017a). We assume that GW localization volumes are spherical, and calculate the radius max gw of the biggest volume depicted in Figure 1. We then populate with AGN a sphere of radius = L (0.2) + max gw , where L (0.2) is the luminosity distance corresponding to = 0.2. The centre of this sphere corresponds to the position of the interferometric network we simulate the detections of. All the AGN are treated as point sources and their distribution is uniform in comoving volume. We then consider a set of gw GW detections and draw for each of them a value of V90 from the relevant distribution in Figure 1. We denote with the localization volume associated to the − ℎ detection. We require that the centre of each has a distance from the interferometric network smaller than L (0.2). A fraction eff agn = 0.9 agn of the centres of the localization volumes are set in order to correspond to the position of an AGN. We here use eff agn instead of agn to take into account the fact that we are here dealing with 90% CL localization volumes, and therefore we expect only the 90% of the origins of the simulated GWs to be actually located in such volumes.
We then count the number of AGN in each localization volume . Equation 1 ensures that each GW localization volume is entirely contained in our simulated Universe.
For every set of gw simulated GW detections, we then calculate where L (0) and L ( agn ) are the likelihood functions of the noconnection hypothesis and of the agn -correlation hypothesis, respectively. These likelihood functions are constructed assuming a Poissonian distribution for . See Bartos et al. (2017a) for more details. Every simulation is therefore associated to a value of that depends on agn , gw , agn , the value of V90 of each simulated GW event, and the number of AGN within such volume.
We expect to be positive in simulations in which gw eff agn localization volumes' centres correspond to an AGN. We refer to simulations that satisfy this requirement as signal realizations, and to the value of obtained from each of them as s . Likewise, we call b every value of that is obtained from a background realization. These realizations are simulations in which the centres of the localization volumes are randomly distributed, uniformly in comoving volume. We, therefore, expect b to be negative.
We perform 3,000 signal realizations and the same amount of background realizations, for each set of values of gw , agn and agn . Once a value for agn and for agn has been set, an increase in gw leads to a greater separation between the distribution of s and the distribution of b .
The target degree of significance in the rejection of the noconnection hypothesis is reached when the median value of the distribution of s corresponds to a p-value lower than 0.00135 when compared to the b distribution.
To evaluate 3 gw for a specific value of agn , we calculate 30 pvalues, keeping such parameter fixed together with gw . We repeat these calculations for multiple values of gw , and then fit the trend of the average p-value for a given gw as a function of gw itself. Such trend is well fitted by a decreasing exponential function for every value of agn we investigated. Once the parameters of these fits are known, we invert the fit function and calculate the number of detections corresponding to a p-value of 0.00135. We repeat the same analysis for 6 different values of agn between 0.5 and 1.

Minimum number of GW detections with fixed agn
In this section, we present the results obtained keeping the AGN number density parameter fixed to agn = 10 −4.75 Mpc −3 . The trend of 3 gw as a function of agn is shown in Figure 2. The error bars correspond to the standard deviation of 1, 000 values of 3 gw calculated for each of the 6 values of agn we test. The results for O3 and O4 are represented by the blue squares and the green dots, respectively.
The trend of 3 gw as a function of agn is fitted with the same functional form used in Bartos et al. (2017a), which is the following: The best fit values we obtain in the case of the O3 simulated events are = 35.8±1.2 and = 1.73±0.08, while for O4 simulated events we obtain = 20.7 ± 0.7 and = 1.73 ± 0.11. We perform the same analysis for O3 with a lower (2 ) significance threshold for the rejection of the no-connection hypothesis. In this case, the best-fit values for the fit are = 17.7 ± 0.5 and = 1.57 ± 0.12.

Significance of the no-connection hypothesis rejection as a function of agn and agn
During the third observing run of the LVK Collaboration, 13 detected BBH mergers have an expectation value of redshift lower than 0.2. As we can infer from the results presented so far, with this low number of "closeby" events it is not possible to reject with a 2 significance the no-connection hypothesis for any value of agn , assuming agn = 10 3 significance the no GW-AGN connection hypothesis as a function of the fraction of GW originated from an AGN. The error bars represent the standard deviation over 1, 000 realizations of 3 gw obtained for each tested value of agn . The results for the third and the fourth observing run of LVK Collaboration interferometers are represented by the blue squares and the green dots, respectively. The data points have been fitted with the following function: 3 agn = − agn . The best-fit values for the O3 scenario are = 35.8± 1.2 and = 1.73 ± 0.08, while for the O4 scenario they are = 20.7 ± 0.7 and = 1.73 ± 0.11. The best-fit function for O3 (O4) is represented by the blue (green) dashed line. becomes more significant, and a lower number of detection is needed to rule out the no-connection hypothesis.
Hence, we perform the same analysis as above but keeping gw fixed at the value of 13, and varying both agn and agn . For each point in this 2D parameter space, we determine the p-value associated to the median of the distribution of when compared to the distribution of .
The results of such analysis are shown in Figure 3. The white dashed (solid) line divide the parameter space into two decision regions, corresponding to parameter choices for agn and agn whose associated p-values are lower or higher than 0.00135 (0.02275), i.e. a significance higher or lower than 3 (2 ), respectively. The noconnection hypothesis can be rejected accordingly in the respective regions.
For example, with 13 GW detections and assuming a number density of AGN agn = 10 −7.50 Mpc −3 , we can, in principle, reject the no-connection hypothesis with a 3 (2 ) significance if agn ≥ 0.40 ( agn ≥ 0.25). Such a low number density corresponds, in the local Universe, to AGN with bolometric luminosities 10 46.2 erg s −1 (Hopkins et al. 2007), or with central SMBHs with masses 10 8.5 M (Greene & Ho 2007).

DISCUSSION AND CONCLUSION
We perform a statistical investigation based on the method presented in Bartos et al. (2017a)    (2) in a simulation in which a fraction agn of GWs come from an AGN. On the other hand, every value of b comes from a simulation in which no GW event is originated in an AGN. On the right-hand axis, we report the logarithm of the minimum bolometric luminosity min [erg s −1 ] that has to be considered in the integration of the AGN L F at = 0.1 (Hopkins et al. 2007) to obtain the value of log( agn [Mpc −3 ]) indicated on the left-hand side of the grid. and GW localization volumes, how many GW detections are needed to reject the no GW-AGN connection hypothesis. We find that the 13 O3 GW detections with expected ≤ 0.2 are not enough to reject the no-connection hypothesis with either 3 or 2 significance. This result is obtained considering AGN with a number density agn = 10 −4.75 Mpc −3 . Nonetheless, Figure 3 shows that with the same number of detections, it is possible to reject the no-connection hypothesis for specific values of the AGN number density and of the fraction of GW events that originated inside an AGN. More precisely, the lower the AGN number density (i.e. the higher the luminosity of the considered AGN, or the higher the mass of the central SMBH), the more likely it is to reject such a hypothesis with a given significance threshold. As far as O4 is concerned, the green line in Figure 2 shows that at least 21 detections associated with redshift ≤ 0.2 will be needed to be able to reject the no-connection hypothesis between BBH mergers and AGN with agn ≥ 10 −4.75 Mpc −3 . The number of expected detections of BBH mergers during O4 is 79 +89 −44 (Abbott et al. 2020). In our simulations of O4 detections, roughly 18.25% of the detected events correspond to ≤ 0.2. Our estimate is therefore that during O4, 14 +17 −8 BBH mergers will be associated with ≤ 0.2. As shown in Figure 2, 30 O4 closeby BBH detections would be enough to test values of agn higher than ≈ 80%, using AGN of any luminosity. The same degree of GW-AGN connection could be tested using a lower number of O4 detections in combination with the 13 closeby O3 detections.
We restrict our analysis to GW events with an expectation value for the redshift of ≤ 0.2 for two different reasons. First, far GW events are typically associated with much larger localization volumes than the ones associated with closer events. The inclusion of large GW localization volumes in our algorithm makes it too computationally demanding. The second reason is that for very luminous AGN in the local Universe, we expect to have high values of completeness in real AGN catalogues. These high values are needed in order to produce reliable results when applying the method described in this work to real, observed GW events and AGN. The incompleteness in observed AGN catalogues can be nonetheless taken into account with an appropriate rescaling of agn (Bartos et al. 2017a).
The main assumptions we made in this work were: considering spherical GW localization volumes, and neglecting redshift evolution for AGN and GW events as well as AGN clustering. We expect these assumptions not to remarkably impact on our final results. The BBH merger rate, the AGN number density, and the expected AGNassisted merger rate do not significantly vary within the redshift range we consider (Hopkins et al. 2007;Yang et al. 2020;The LIGO Scientific Collaboration et al. 2021c). Taking into consideration the real shape of GW events localization volumes and the clustering of AGN in the local Universe is important when performing a maximum likelihood estimation to find which value of agn best represents real observations. Such estimation is not the aim of this work but is currently being implemented in ongoing projects, in which the exploitation of realistic AGN catalogues and GW skymaps is required.
Our results motivate more in-depth statistical investigations of a possible connection between GW events and AGN exploiting actual data. An quantitative assessment of agn would allow us to put constrain on the BBH merger rate per AGN in the local Universe, that in turn can be used to both deepen our physical understanding of AGN disks and inform theoretical models of BBH mergers from the AGN formation channel. Finally, extrapolating these findings at higher redshift will enable predictions for GW events detectable by the third-generation GW interferometers.