Gravitational Wave Source Populations: Disentangling an AGN Component

We construct a one-parameter model for binary black hole (BBH) formation and merger within an AGN disk, parameterized by the maximum mass ${m}_{\max }$ of the natal black hole (BH) distribution. Specifically, we adopt a seed BH mass distribution which follows the Salpeter mass function with index 2.35, ${dN}/{dm}\propto {m}^{-2.35}$ with given m_max. We considered different m_max cases to probe the multiple AGN models (pushing m_max in the nucleus higher) along with probing the history of the nuclear cluster. For neutron stars (NSs) we assume a normal distribution m/M_⊙ ∼ N(1.49, 0.19). The BHs and NSs are assumed to orbit a supermassive black hole in an AGN, migrating into the disk and inward from their natal locations. Close to the AGN, these objects undergo multiple encounters, facilitating binary formation and merger. Other AGN parameters are fiducial values that are expected to be typical, while there are large uncertainties. The range of possible values in the AGN models parameter space is discussed in McKernan et al. (2018).

Following Bartos et al. (2017), we adopted a geometrically thin, optically thick, radioactively efficient, steady-state accretion disk expected in AGNs. We used a viscosity parameter α = 0.1, radioactive efficiency = 0.1, fiducial supermassive BH mass M_• = 10⁶ M_⊙ and accretion rate $0.1{\dot{M}}_{\mathrm{Edd}}$, where ${\dot{M}}_{\mathrm{Edd}}$ is the Eddington accretion rate. Using Yang et al. (2020b) and Tagawa et al. (2020b, 2020a), we have computed the expected mass and spin distributions of binary mergers in AGNs.

Figure 1 shows the binary black hole merger intrinsic parameter distribution for the AGN model with different initial mass limits (${m}_{\max }=[15{M}_{\odot },35{M}_{\odot },50{M}_{\odot },75{M}_{\odot }]$). As the natal BH mass upper limit m_max increases, more-massive binary components and total masses are allowed. At the same time, for high m_max, asymmetric systems are more frequent. By contrast, we have observed that the spin distribution properties are largely independent of our choice for the m_max limit.

Zoom In Zoom Out Reset image size

To allow for binary black holes that have a non-AGN origin, we follow previous work and introduce a few-parameters mixture model family. As illustrated in Figure 2, for the non-AGN component, we allow binary black holes to arise from a mixture of a power law and Gaussian components, as detailed in Table 11 of Abbott et al. (2021a). In the power-law component, the primary is drawn from a pure power-law mass distribution (with some unknown ${m}_{\max ,\mathrm{pl}}\lt 50{M}_{\odot }$ and unknown primary power-law index); the mass ratio is drawn from another power law; and the spins are drawn from an unknown beta distribution. In the Gaussian component, the primary and secondary are drawn from two independent Gaussian distributions with unknown mean and variance, with both means confined a priori to be near 30–40 M_⊙ to be consistent with expectations for PISN supernovae. By including a second component this phenomenological model can allow for non-power-law features and also allow for spin distributions that vary with mass. Each component k has some undetermined overall rate ${{ \mathcal R }}_{k}$. The top panel of Figure 2 shows the general non-AGN model.

Zoom In Zoom Out Reset image size

Our overall model is therefore a mixture model, parameterized by the unknown (continuous) AGN merger rate and its (discrete) maximum mass ${m}_{\max }$, along with all parameters of the non-AGN mixture model. The overall merger rate density dN/dVdXdt can therefore be expressed as a sum

$$\begin{eqnarray*}&&\displaystyle \frac{{dN}}{{dVdXdt}}={{ \mathcal R }}_{\mathrm{agn}}{p}_{{agn}}(X| {m}_{\max })+{{ \mathcal R }}_{g}{p}_{g}(X| {{\rm{\Lambda }}}_{g})+{{ \mathcal R }}_{\mathrm{pl}}{p}_{\mathrm{pl}}(X| {{\rm{\Lambda }}}_{\mathrm{pl}}),\end{eqnarray*}$$

where X are binary parameters and where p_q , Λ_q are the model distributions and parameters for the qth component (AGN, Gaussian, and power-law respectively).

To systematically assess how well the distinctive features in AGN formation scenarios can be disentangled from this large model family, we will perform a sequence of calculations with increasing model complexity, as shown in the panels of Figure 2. Specifically, we consider without any non-AGN component; the power-law and Gaussian (PL+G) model only; the power-law and AGN model (PL+AGN), without any Gaussian component; and finally the most general model with all three components.

We describe and demonstrate a flexible parametric method to infer the event rate as a function of compact binary parameters, accounting for Poisson error and selection biases. In (Abbott et al. 2021b), analyzed the multispin model, which is the joint mass-spin model for binary black holes. Independent analyses have shown that there is a feature in the BBH mass spectrum around 30–40 M_⊙, which is modeled as a Gaussian peak on top of a power-law continuum. It is an empirical way of modeling extra features, but here we tried to understand it feature with AGN models.

In this section we review the population inference Wysocki et al. (2019) used for multiformation channel contributions. Binaries with intrinsic parameters x would merge at a rate ${dN}/{{dV}}_{c}\,{dtdx}={ \mathcal R }\,p(x)$, where N is the number of detections, V_c is the comoving volume, ${ \mathcal R }$ is the spacetime-independent rate of binary coalescence per unit comoving volume and p(x) is the probability of x from detected binaries. The binaries intrinsic parameters includes mass m_i and spins S _i , where i = 1, 2.

The likelihood of the astrophysical BBH population at a given merger rate ${ \mathcal R }$ (Loredo 2004; Mandel et al. 2019; Thrane & Talbot 2019) and given binary intrinsic parameters X ≡ (m₁, m₂, χ ₁, χ ₂), where ${{\boldsymbol{\chi }}}_{i}={{\boldsymbol{S}}}_{i}/{m}_{i}^{2}$, given the data for N detections ${ \mathcal D }=({d}_{1},\ldots ,{d}_{N})$. The likelihood is given by

$$\begin{eqnarray}&&{ \mathcal L }({ \mathcal R },X)\equiv p({ \mathcal D }| { \mathcal R },X),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{ \mathcal L }({ \mathcal R },X)\propto {e}^{-\mu ({ \mathcal R },X)}\,\displaystyle \prod _{n=1}^{N}\int {dx}\,{{\ell }}_{n}(x)\,{ \mathcal R }\,p(x,X),\end{eqnarray} \tag{ 2 }$$

where $\mu ({ \mathcal R },X)$ is the expected number of detections under a given population parameterization X. Using the Bayes theorem, one may obtain a posterior distribution on ${ \mathcal R }$ and X after assuming some prior $p({ \mathcal R },X)$. For computational efficiency and to enable direct comparison with discrete formation models, we use the Gaussian likelihood approximation technique introduced in Delfavero et al. (2022) to characterize each BBH observation's likelihood ℓ_n (x).

Using this formalism, we estimate what fraction of binary black holes are generated via the AGN channel using detected gravitational wave merger information. To do this study we have upgraded current parametric methods with a mixture model feature. Here we have a freedom to do the analysis with a number of models, which has different astrophysical binary distributions.

We have estimated the astrophysical black hole merger rate for a given astrophysical model with a given number of gravitational wave detections. For this study, we have used the parameter estimation samples obtained by the LIGO–Virgo–KAGRA Collaboration (Abbott et al. 2019, 2021c, 2021d, 2021g) and those available in the Gravitational Wave Open Science Center ⁶  with the mixed model. Here we follow the same selection criteria as Abbott et al. (2021f); we have considered events with a false alarm rate of <0.25 yr⁻¹.

Table 1 shows the merger rates inferred for the joint power-law (PL)-only, G-only, AGN-only (with different m_max), and PL+AGN model (with different m_max) and PL+G+AGN models (with different m_max). We have estimated the merger rate results derived using each individual model component as well as combined. We have observed that the inferred merger rate from the AGN-only models with different choices for ${m}_{\max }$ produces largely consistent results peaking near 50 Gpc⁻³ yr⁻¹. For the smallest maximum BH natal mass ${m}_{\max }=15{M}_{\odot }$, the inferred single-component AGN merger rate peaks around 80 Gpc⁻³ yr⁻¹. For all other models, they peak around the same ${ \mathcal R }$. Similarly, we have inferred the merger rate distribution for the single-component Gaussian, and power-law components. As expected, the merger rate from the power-law component dominates overall, as it incorporates and describes many frequent mergers of the lowest-mass binary black holes. In the case of PL+AGN, the inferred AGN and PL components have a well-determined merger rate, with a median AGN merger rate ≃8/,Gpc⁻³ yr⁻¹ and a PL merger rate ≃30/,Gpc⁻³ yr⁻¹. Changing the maximum natal BH mass m_max has a very mild impact on the inferred AGN merger rate, and almost no effect on the inferred PL merger rate.

Table 1. The Astrophysical Rates from PL+AGN and PL+AGN+G Models

Analysis	Models	${m}_{\max }=15$	${m}_{\max }=35$	${m}_{\max }=50$	${m}_{\max }=75$
PL only	${71.9}_{-20.4}^{+19.9}$	⋯	⋯	⋯	⋯
G only	${19.1}_{-2.3}^{+2.5}$	⋯	⋯	⋯	⋯
AGN only		${84.7}_{-18.5}^{+19.5}$	${49.8}_{-5.5}^{+6.1}$	${52.9}_{-6.1}^{+7.8}$	${53.2}_{-6.2}^{+7.7}$
PL+AGN	PL	${29.3}_{-7.2}^{+9.9}$	${28.8}_{-6.8}^{+11.3}$	${25.7}_{-6.3}^{+7.6}$	${26.8}_{-6.2}^{+9.5}$
	AGN	${8.7}_{-4.5}^{+6.3}$	${8.3}_{-3.7}^{+8.3}$	${12.7}_{-4.1}^{+6.5}$	${11.9}_{-3.2}^{+4.6}$
PL+AGN+G	PL	${21.3}_{-5.2}^{+7.3}$	${23.5}_{-5.5}^{+6.6}$	${21.6}_{-5.0}^{+6.7}$	${23.3}_{-5.0}^{+7.1}$
	AGN	${6.7}_{-4.0}^{+5.6}$	${6.2}_{-3.5}^{+5.0}$	${12.7}_{-4.9}^{+5.0}$	${12.9}_{-4.3}^{+4.3}$
	G	${10.2}_{-2.3}^{+2.2}$	${9.3}_{-3.7}^{+2.6}$	${5.9}_{-1.9}^{+2.9}$	${5.6}_{-1.5}^{+2.2}$

Note. Each row corresponds to each analysis and each column corresponds to different AGN models with different initial mass limits.

Download table as:  ASCIITypeset image

In the case of PL+G+AGN analyses, the inferred AGN, G, and PL components have a well-determined merger rate, with a median AGN merger rate ≃7/,Gpc⁻³ yr⁻¹, a G merger rate ≃4/,Gpc⁻³ yr⁻¹, and a PL merger rate ≃30/,Gpc⁻³ yr⁻¹. As we have seen in the PL+AGN study, we have not seen any major effect on PL or G merger rate when we change the AGN model. Our inferred merger rates deduced from the multicomponent model are consistent with inferences performed using single-component models alone, suggesting that inference isolates only the contribution from each component.

To better appreciate how well our inference directly projects out the relative contribution from each component, Figures 3 and 4 show our inferred merger rate versus mass for PL+AGN and PL+G+AGN models analyses respectively. In each plot we have shown each model component in mass space.

Zoom In Zoom Out Reset image size

As we expected, the low-mass region is highly contributed by the PL model compared to other models. We have observed a peak in PL model distribution around 7M_⊙ for both PL+AGN as well as PL+G+AGN analyses. The peak is prominent for the PL+G+AGN study compared to PL+AGN. The high-mass region is represented by only AGN, as we expected to see. Note that the AGN model not only contributes in the high-mass region, it also contributes in full mass space, as shown in Figure 3.

Zoom In Zoom Out Reset image size

Similarly here we show the inferred merger rate versus mass ratio for PL+AGN and PL+G+AGN models analyses.

Figures 3 and 4  show our inferred merger rate versus mass ratio for PL+AGN and PL+G+AGN models analyses respectively. In each plot we have shown each model component in mass ratio space. As we expected, the AGN model contributes to the low-mass ratio region for PL+AGN analysis and the AGN and G models do so for PL+G+AGN analysis. For high q, the dominate contribution was from PL model, which is consistent with detected events.

As we discussed before, the contribution of a power-law model to the overall merger rate does not change substantially if we include or omit other model components like AGN or G. Among the models we consider, this quasi-universality is expected: the PL model most effectively reproduces the merger rate versus mass for the lowest-mass and most frequently merging binary black holes. While the overall merger rate from this component is stable to our choice of the mixture, the model parameters recovered for PL depend strongly on which other confounding contributions are also present, as suggested by Figures 3 and 4. While the PL mass ratio distribution does not depend strongly on including or omitting AGN or G, the power-law slope α and minimum mass m_min do change substantially. The estimated α median value with 65% credible intervals from different analysis as ${1.6}_{-0.2}^{+0.2}$, ${6.7}_{-2.4}^{+3.1}$,${8.4}_{-3.4}^{+2.3}$, and ${1.8}_{-0.5}^{+0.5}$ for PL-only, PL+G, PL+AGN (m_min = 50), and PL+G+AGN (m_min = 50) respectively. Similarly, m_min estimation are ${2.5}_{-0.3}^{+0.3}$, ${8.4}_{-0.3}^{+0.2}$, ${8.5}_{-0.4}^{+0.2}$, and ${6.7}_{-1.3}^{+1.6}$ for PL-only, PL+G, PL+AGN (m_min = 50), and PL+G+AGN (m_min = 50) respectively. For example, the α estimation suggests that while a pure power-law model favors ${m}_{\min }$ close to the lower limit our priors allow, incorporating other components causes the power-law component's minimum mass to favor larger masses. With the pertinent mass range for the power law changing substantially via different ${m}_{\min }$, unsurprisingly. The α estimation has a wide range of inferred power-law exponents, as the PL may dominate only an extremely narrow range of masses; see Figure 3.