Measuring the properties of active galactic nuclei disks with gravitational waves

Active galactic nuclei (AGN) are promising environments for the assembly of merging binary black hole (BBH) systems. Interest in AGNs as nurseries for merging BBH is rising following the detection of gravitational waves from a BBH system from the purported pair-instability mass gap, most notably, GW190521. Active galactic nuclei have also been invoked to explain the formation of the high-mass-ratio system, GW190814. We draw on simulations of BBH systems in AGN to propose a phenomenological model for the distribution of black hole spins of merging binaries in AGN disks. The model incorporates distinct features that make the AGN channel potentially distinguishable from other channels, such as assembly in the field and in globular clusters. The model parameters can be mapped heuristically to the age and density of AGN disks. We estimate the extent to which different populations of mergers in AGNs can be distinguished. If most merging black holes are assembled in AGNs, future gravitational-wave observations may provide insights into the dynamics of AGN disks.


INTRODUCTION
Gravitational-waves from the mergers of binary black hole (BBH) systems have recently transformed astronomy. However, the astrophysical origins of these events are still uncertain. There are two main proposed astrophysical pathways to the mergers: (i) isolated binary evolution via mass transfer, including a common envelope phase in galactic fields; (ii) dynamical formation in dense environments. Each pathway is associated with different distributions of black-hole spin (Mandel & O'Shaughnessy 2010;Stevenson et al. 2017;Fishbach et al. 2017;Talbot & Thrane 2017;Wysocki et al. 2019) and binary eccentricity (Rodriguez et al. 2018;Samsing 2018;Lower et al. 2018;Romero-Shaw et al. 2019;Romero-Shaw et al. 2021;Zevin et al. 2021;Gayathri et al. 2020). Measuring BBH spins and eccentricity with gravitational waves can therefore be used to determine how and where BBH are assembled (Abbott et al. 2019;Abbott et al. 2021).
Active Galactic nuclei (AGNs) are expected to contain a dense population of stars and stellar remnants, avi.vajpeyi@monash.edu such as stellar-origin black holes (BHs) (Morris 1993;Miralda-Escudé & Gould 2000;Hailey et al. 2018;Generozov et al. 2018). Binary black hole systems can form via close encounters in this dynamically "hot" environment but are often rapidly "ionized" via tertiary encounters (Antonini & Rasio 2016;Fragione et al. 2019). The dense nuclear population and AGN gas disks (when present) can interact, resulting in an embedded population of stars and BH within the disk. These embedded objects can weakly perturb the surface-density profile of the gas disk, resulting in gas torques within AGN disk that allow for Type I (non gap-opening) migration of the embedded objects (McKernan et al. 2012). Differential migration rates of the objects encourage binary formation, leading to compact binary mergers detectable with LIGO-Virgo (McKernan et al. 2014; Bartos et al. 2017;Stone et al. 2017;Wang et al. 2021b;LIGO Scientific Collaboration et al. 2015;Acernese et al. 2015). Kicked merger products are generally retained by the deep potential well, allowing for hierarchical BBH mergers. If most mergers observed by LIGO-Virgo were assembled in an AGN disk, it may be possible to reverse-engineer conditions beneath the AGN photosphere Wang et al. 2021b).
The gaseous disk in AGNs likely serve to align (to varying degrees) both black-hole spin vectors 1 χ 1,2 and the binaries' orbital angular momentum vectors Ldepending on the density and age of the disk (Bogdanović et al. 2007). On the other hand, tertiary encounters with binaries in the disk tend to misalign χ 1,2 relative to L-depending on the timescales of encounters (Liu & Lai 2017;Tagawa et al. 2020a). The competing effects of the gaseous disk and dynamical encounters on BBHs in AGNs determine the distribution of BBH spin orientations. On the other hand, binaries born in the field have χ 1,2 nearly aligned to L with a small spread due to supernova kicks (Kalogera 2000;Mandel & O'Shaughnessy 2010;Dominik et al. 2013;Giacobbo et al. 2017;Eldridge et al. 2017;Olejak et al. 2020). Finally, dynamically assembled binaries (e.g., in globular clusters) exhibit no correlation between the L and χ 1,2 . In this paper, we propose a phenomenological model for the distribution of black-hole spins in AGNs to capture the salient features predicted from theoretical modeling.
The remainder of this paper is organised as follows. In Section 2, we review the spin-orientations of BBH in AGNs. Section 3.1 presents a phenomenological model describing AGN BBH spins-orientations. In Section 4 we present the results of a simulated study, demonstrating the model's effectiveness.

SPIN PROPERTIES OF BINARY BLACK HOLES AT FORMATION
Modeling AGNs is challenging due to the interplay between gas dynamics, scattering binaries, and feedback from the central supermassive BH. Commonly used disk models span wide ranges of disk density and geometry (see, e.g., Sirko & Goodman 2003;Thompson et al. 2005), and provide broad estimates for merger rates in AGNs, (e.g., McKernan et al. 2018;Gröbner et al. 2020;Tagawa et al. 2020a). However, we have qualitative predictions for the spin-orientation population properties of merging BBH systems using Monte-Carlo and N -body simulations to identify key features. Some of these predictions are tabulated in Table 1.
In this section, we discuss three predictions for the spin distribution of merging BBH at the time the binary is formed 2 . The spin vectors subsequently evolve through general relativistic precession of the orbital plane. Nonetheless, the orientation of the spin vectors at the merger contains information about the orientation at formation. We describe the phenomenology of the spin vectors at the formation and then discuss later the resulting phenomenology at the time of the merger.
(i) Gas accretion torques BH spins to align with the disk. The BH embedded in the AGN disk early in the disk lifetime should have isotropically distributed spins at formation (Mapelli & Gualandris 2016;Tagawa et al. 2020a). As the BH of mass m migrate in the disk, they accrete ∆m disk-gas mass, resulting in a torque pointing into the plane of the AGN disk. The magnitude of the torque on BH spins depends on ∆m (Bogdanović et al. 2007) and can be summarised as follows.
Dense gas disks can torque initially randomly oriented BH spins into alignment with the disk angular momentum in < 5 Myr (e.g. McKernan et al. 2020), assuming an Eddington-limited accretion rate. Fabj et al. (2020) find that a critical density of ρ > 10 −11 gcm −3 is required to capture many orbiters into the disk over a modest AGN lifetime. So, here we assume that gas disks with densities ρ < 10 −11 gcm −3 need to live > 5 Myr both to capture enough BH over their lifetimes to make a significant contribution to the BBH merger rate and to torque BH spins into alignment with the disk.
In sufficiently long-lived, dense disks (Table 1a), fully embedded BH accrete more than 1 − 10% of their initial mass ∆m 0.01 m i . The resultant torque from gas accretion onto the embedded BH reorients the BH spin vector χ to align with the angular momentum vector for the AGN disk J AGN , (e.g., Bogdanović et al. 2007).
Note that the details of accretion onto objects embedded in AGN disks is subject to much uncertainty. Feedback, turbulence, and interactions can alter gas flow dynamics inside the BBH Hill sphere, possibly inhibiting a high accretion rate on component BH, limiting the average torque magnitude (Hankla et al. 2020). Lower accretion rates can result in longer timescales (τ AGN 5 Myr) for BH being torqued into alignment with the disk. Table 1. The phenomenology of black holes spin for binaries merging in AGNs. The two rows correspond to different values of the population parameter σ1 defined in Section 3.1, which controls the effective density of the AGN. The two columns correspond to different values of the population parameter σ12 defined in Section 3.1, which controls the effective age of the AGN. Each cell is divided into two. The cell's left-hand side describes the AGN's properties while the right-hand side describes the distribution of black hole spins at the time of formation.
(ii) Gas torques dampen BBH orbital angular momentum. Binary black hole systems in AGN disks experience Lindblad and co-rotating gas torques as they migrate through the disk (Tagawa et al. 2020b). Co-rotation torques dampen the binary's eccentricity and drive L into alignment with J AGN on a characteristic timescale (Tanaka et al. 2002) where M SMBH is the SMBH mass, m b is the BBH system's total mass, h = H/a the disk aspect ratio, Σ is the disk surface density and Ω is the Keplerian orbital frequency. In addition to aligning L with J AGN , dynamical gas friction can promote binary hardening-the process of losing orbital energy and tightening the orbit (Baruteau et al. 2011). If t damp is smaller than the lifetime of the AGN (typically in long-lived AGN τ AGN 5 Myr, Table 1ac) the BBH orbital angular momentum vector L will align with J AGN by the time the binary merges.
Multiple migrators can quickly interact with each other, potentially leading to complex or chaotic dynamical encounters moderated by the disk gas (Wang et al. 2021a). Tertiary encounters of binaries with compact objects in the disk or the spherical nuclear population component can harden or soften BBH systems (Leigh et al. 2018;Yang et al. 2019a), increase BBH orbital eccentricity (Samsing et al. 2020), and alter the orbital angular momentum of the BBH (Tagawa et al. 2020b,a). Close encounters with a tertiary object on a diskcrossing orbit can perturb the orbital angular momentum of the BBH on the timescale for the encounter (t enc Leigh et al. (2018)), which depends on the density of the nuclear star cluster (ρ NSC ), BBH location in the disk, and the efficiency of disk capture (a function of ρ, τ AGN ). Small values of t enc (Table 1cd) lead to more binaries with L misaligned with J AGN at merger.
The relation between the dampening and dynamical encounter timescales may provide further details about the AGN. For example, t damp > t enc could occur if (I) AGN disks are not long-lived τ AGN ∼ 0.5 − 5 Myr, and the spherical population is not efficiently captured by the disk, resulting in lots of dynamic interactions (McK-ernan et al. 2018;Tagawa et al. 2020b), (II) BBH are positioned in a short-lived inner disk where the encounter rate with the spherical component is significantly higher (a is small, Leigh et al. (2018)), (III) ρ NSC is large (e.g., the NSC is cuspy not cored McKernan et al. (2018); Tagawa et al. (2020b)), or (IV) The fraction of BBH systems hardened via gas torques gas 3 is not particularly efficient relative to dynamical hardening (Stone et al. 2017).
In summary, our assumptions are that: 1) BH embedded in AGN disks accrete at the Eddington rate and are thereby gradually torqued into alignment with the disk, 2) BH from the spheroid are captured by the disk over its lifetime (preferentially heavier BH at lower inclination angles), such that heavier BH tend to spend more time embedded in the AGN disk and therefore experience longer periods of torquing. Torquing into alignment depends on the gas mass accreted (i.e., a combination of the disk gas density and the disk lifetime for a given accretion rate) and how long a BH has been in the disk. For example, a short-lived but dense gas disk can more rapidly torque embedded BH into alignment than a less dense gas disk. However, BH that are captured by such a disk will spend less time embedded and will experience a shorter duration torque towards alignment.
From such assumptions, old (i.e., long-lived), dense AGNs (Table 1a) produce a population of BBH systems with χ 1 preferentially aligned with the orbital angular momentum and χ 2 preferentially aligned with χ 1 . Midaged, dense AGNs (Table 1b) produce BBH systems with χ 1 preferentially aligned with the orbital angular momentum, but χ 2 is not correlated with χ 1 . This is because we assume M 1 spends more time embedded in the AGN disk on average. Old, dilute AGNs (Table 1c) produce BBH systems where χ 1 is not preferentially aligned with the orbital angular momentum. However, χ 2 is preferentially aligned with χ 1 . Young, dilute AGNs (Table 1d) produce BBH systems where χ 1 is uncorrelated with the orbital angular momentum and χ 2 is uncorrelated with χ 1 .

Model Description
In this section, we construct a model for the spin orientation of black holes (at the time of formation) in merging binaries residing in the AGN disk. It will be useful to employ two coordinate systems. Vectors with no prime are measured with respect to the orbital an-gular momentum vector such that while primed vectors are measured with respect to the primary spin vector such that Hence, χ 1 is given by in the unprimed coordinate system and in the primed coordinate system. A schematic diagram depicting χ 1 , χ 2 , L and the angles θ 1 , θ 12 can be seen in Figure 1. Figure 1. Schematic diagram of a binary black hole system. The blue, orange, and purple arrows correspond to χ 1, χ 2 and L . The angle between ( χ 1, L ) is θ1, while the angle between ( χ 1, χ2 ) is θ12.
Here, θ 1 is the zenith angle between L at formation andχ 1 , and φ 1 is the azimuthal angle measured fromx, aboutẑ at formation. The secondary spin vector is where θ 12 is the zenith angle between the χ 1 and χ 2 , and ζ 12 is the azimuthal angle measured fromx , about χ 1 .
The hyper-parameters determine the shape of the distribution.
This distribution allows for preferred alignment between χ 1 and L with a free parameter σ 1 controlling the typical misalignment angle. For small values of σ 1 , thê χ 1 distribution tends to be nearly aligned with the L.
As σ 1 becomes large, the distribution becomes uniform, soχ 1 is uncorrelated with L. We assume the primary azimuthal spin angle φ 1 is drawn from a uniform distribution denoted U . We assume that the secondary spin vector is preferentially aligned to the primary spin vector 4 at formation such that π(cos θ 12 |σ 12 ) = N t (cos θ 12 |σ 12 ).
Small values of σ 12 imply that χ 1 and χ 2 tend to point in nearly the same direction. As σ 12 becomes large, the directions of χ 1 and χ 2 become uncorrelated. We assume that ζ 12 is drawn from a uniform distribution.
3.2. The evolution of spin vectors with time 4 We take the phrase 'preferentially aligned spin' to mean that the directions of χ 1 and χ 2 vectors are correlated so that they point more nearly in the same direction than two random vectors.
As BBH systems evolve, coupling of the component black hole spin vectors and the orbital angular momentum vector results in Lense-Thirring precession, causing the spin and angular momenta to precess about the system's total angular momentum vector (Mashhoon et al. 1984). Thus, the angles θ 1 , θ 12 evolve over time. While the distributions of θ 1 , θ 12 modeled in Eq. 10 describe the binary properties at formation (approximately when the binary is infinitely separated, Johnson-McDaniel et al. (2021)), these distributions are different by the time these binaries enter the observing band of audioband gravitational-wave detectors (20-20,000 Hz).
Fortunately, the information present in the distribution of θ 1 , θ 12 at formation is encoded in the distribution of the "effective inspiral spin parameter" χ eff (Damour 2001) and the effective precession parameter χ p (Schmidt et al. 2012), which are approximate constants of motion (Hannam et al. 2014;Gerosa et al. 2021). The parameter χ eff measures the spin components aligned with the orbital angular momentum while χ p measures the spin components in the orbital plane. Figure 3 shows joint distributions of χ eff , χ p , each representing three different AGN populations, each with different values of σ 1 , σ 12 . The dashed-orange distribution is created using σ 1 = 0.1 and σ 12 = 10 (a mid-aged, dense AGN). The solid-purple distribution displays the distribution for σ 1 = 0.1 and σ 12 = 0.1 (an old, dense AGN). Finally, the solid-green contours is the distribution with σ 1 = 10 and σ 12 = 10 (a young, dilute AGN). The distinguishability of these distributions illustrates how the AGN properties may be imprinted on the distribution of quantities measured by LIGO-Virgo. To recast our model in terms of quantities that are measured by LIGO-Virgo, it is necessary to compute π(χ eff , χ p |σ 1 , σ 12 ). (11) In principle, an expression for this distribution may be obtained through a series of convolutional integrals, which most likely have to be evaluated numerically. An alternative to numerical integration is estimating the probability density distribution at fixed values of σ 1 , σ 12 using histograms or kernel density estimators. These estimates can be used to interpolate the probability density for arbitrary values of σ 1 , σ 12 , for example, using a machine learning algorithm; see Hernandez Vivanco et al. (2020). Alternatively, a machine learning algorithm can likely be trained to reproduce the results of numerical integrals for different values of σ 1 , σ 12 . However, as this work focuses on estimating how well LIGO-Virgo will be able to measure σ 1 , σ 12 , creating a machine-learning representation of Eq. 11 goes beyond our present scope. Instead, in the demonstration that follows in Section 4, we pretend that the values of θ 1 , θ 12 measured by LIGO-Virgo are unchanged since the binary was formed. This is tantamount to assuming that the χ eff , χ p distribution encodes-without loss-all of the information in the distribution of θ 1 , θ 12 . Since there is likely information loss, our results are overly optimistic. Nonetheless, judging by the distinguishability of different populations shown in Fig. 3, we believe this assumption still yields a ballpark estimate.
In Table 1, we summarize the implications of measuring various σ 1 and σ 12 on some parameters describing AGNs. Some of the implications in Tab. 1 are degenerate. For example, if AGN disks are relatively low density on average (ρ < 10 −11 gcm −3 ), then their corresponding average lifetimes τ AGN 5 Myr to torque χ 1,2 into alignment with J AGN . Moreover, the top left quadrant (small σ 1 , σ 12 ≤ 1, Table 1a) is indistinguishable from a population of BBH systems assembled in the field binary (aligned spins). The bottom right quadrant (large σ 1 , σ 12 > 1, Table 1d) is consistent with a dynamical assembly in a dense stellar environment such as a globular cluster (isotropic spins). Thus, depending on the nature of BBH assembly in AGN disks, this framework does not necessarily provide a useful means of testing the AGN hypothesis against competing hypotheses. Rather, it is a means of probing AGN physics assuming the AGN hypothesis is true. Independent evidence, e.g., from electromagnetic counterparts, may be required to establish this premise.

DEMONSTRATION
In this section, we carry out a demonstration analysis to estimate the ability of Advanced LIGO-Virgo to measure AGN physics. We generate two simulated populations of one hundred BBH signals drawn from our population model described by Eq. 10. The BBH systems are uniformly distributed in comoving volume. We employ the mass model from Talbot & Thrane (2018) with parameters consistent with results from Abbott et al. (2021). For the sake of simplicity, we assume fixed values for the mass-ratio q = 1 and dimensionless spin magnitude χ 1 = 0.6, χ 2 = 0.6. These simplifying assumptions roughly match the mass-ratio and spin magnitudes of the 20% of merging BBH systems with non-zero spin (Roulet et al. 2021;Galaudage et al. 2021). We additionally use an arbitrary polarisation angle φ = 0.1 to simplify analysis further. We draw the black hole spin tilts for each population from: Population A: π AGN ( χ 1 , χ 2 | σ 1 = 0.5, σ 12 = 1), Population B: π AGN ( χ 1 , χ 2 | σ 1 = 3.0, σ 12 = 0.25) at a reference frequency of 20 Hz 5 . As discussed in Section 3.2, we ignore the evolution of the spin vectors from the time of formation to the moment they enter the LIGO-Virgo band, which is equivalent to assuming that the information in the θ 1 , θ 12 distribution is encoded in the distribution of χ eff , χ p without any data loss. This assumption makes our results overly optimistic. Population A corresponds to mergers from AGNs that are old and dense (Table 1a). On the other hand, Population B corresponds to mergers from AGNs that are old and dilute (Table 1c).
We simulate the gravitational-wave signals from the BBH mergers using the waveform approximant IMRPhenomXPHM (Pratten et al. 2020(Pratten et al. , 2021García-Quirós et al. 2020). We add the simulated signals into Gaussian noise coloured to the Advanced LIGO design sensitivity for the Hanford and Livingston detectors (LIGO Scientific Collaboration et al. 2015;Abbott et al. 2018). We ensure that the matched-filter signalto-noise ratio SNR ≥ 8. We perform Bayesian inference with parallel-bilby (Ashton et al. 2019b;Smith et al. 2020;Romero-Shaw et al. 2020;Speagle 2020;Skilling 2004Skilling , 2006 and GWPopulation (Talbot et al. 2019) to recover posterior probability densities for the parameters of the simulated signals using the same waveform approximant IMRPhenomXPHM, also evaluated at a reference frequency of 20 Hz. The prior and population probability distributions we use are documented in Table 2.
In Fig. 4 we plot cumulative distributions of cos θ 1 , cos θ 12 for Population A (old and dense, blue) and Population B (old and dilute, red). The true population distributions are the solid dark curves. The empirically observed distributions from the simulated catalogs are the shaded bands. The darker region of the shaded bands indicates the 90% credibility range of the observed distributions, while the lighter region indicates the 99% credibility range. Each row of the figure displays the observed distributions for different sizes of the populations. As the number of events in each population increases, the observed 99% credibility ranges shrink. Figure 4 demonstrates that the two populations can be visually distinguished once there are O(10) gravitational-wave events. However, given recent results that suggest only ≈ 20% of BBH systems contain a black hole with mea-5 We use 20 Hz for consistency with GWTC-2 (LIGO Scientific Collaboration et al. 2021), a frequency close to the lower-end of the LIGO sensitivity curve (LIGO Scientific Collaboration et al. 2015;Abbott et al. 2018). Note that a reference frequency will not be required when switching to π(χ eff , χp|σ 1 , σ 12 ) from π(cos θ 1 , cos θ 12 |σ 1 , σ 12 ). surable spin (Roulet et al. 2021;Galaudage et al. 2021;Miller et al. 2020)-and given the loss of information moving from θ 1 , θ 12 to χ eff , χ p -it likely that O(50) events are required.

DISCUSSION
This paper introduces a physically motivated phenomenological population model describing the atformation spin-orientations of merging binary black holes assembled in active galactic nuclei. By measuring the distribution of spin orientations with gravitational waves, we may be able to learn about the AGN environment; whether it is old or young, whether it is dilute or dense. We demonstrate that O(10) gravitationalwave events from BBH mergers with SNR ≥ 8 and spin magnitudes χ ∼ 0.6 are required to infer the population parameters describing the shape of the cos θ 1 and cos θ 12 distributions. However, since the majority of BBH mergers include black holes with negligible spin (Roulet et al. 2021;Galaudage et al. 2021;Miller et al. 2020)-and since we have optimistically assumed that the distribution of χ eff , χ p retain all the information of Table 2. Population and Prior distributions for demonstration. We use shorthand to represent distributions: δ delta, U uniform, CM uniform in comoving volume, C cosine, and finally Nt truncated normal (mean of 1). The definitions of the parameters are documented in Romero-Shaw et al. (2020, Table E1). The 'Population' and 'Prior' columns display the spin orientation distributions for the population and prior. The 'Common' displays some parameter distributions common between the population and prior.
To derive this estimate, we assume that all BBH mergers take place in AGNs. This is a reasonable starting point since it is desirable to see if all BBH detections can be understood within a single channel, and it would be somewhat surprising if it turns out that two different channels produce comparable merger rates. However, it may turn out that AGNs provide only some fraction of the observed population. If so, our effort to infer the properties of AGN disks will be complicated by contamination from other channels. One could use the model proposed here as a sub-population in a mixture model, but this goes beyond our present scope.
Our model is cast in terms of the orientation of spin vectors at the time of BBH formation. While the spin vectors subsequently evolve due to precession, the initial orientation of spin vectors is imprinted on the effective spin parameters χ eff , χ p ; see Fig 3. To further develop this model such that it can be used for Bayesian inference, it will be necessary to recast the model in terms of these effective parameters (see Eq. 11). Future work will focus on developing a computationally efficient representation of this distribution, for example, using a machine learning algorithm. With a computationally efficient model, it will be possible to apply the model to current gravitational-wave catalogs, which may now have enough events to begin resolving properties of AGN physics.

ACKNOWLEDGMENTS
We thank Shanika Galaudage and Colm Talbot for their technical help. We thank Ilya Mandel for their helpful advice.
We gratefully acknowledge the Swinburne Supercomputing OzSTAR Facility for computational resources. All analyses (including test and failed analyses) performed for this study used 130K core-hours on OzS-TAR. This would have amounted to a carbon footprint of ∼ 7.6 t CO 2 (Australian Government -Department of the Environment and Energy 2021). However, as OzS-TAR is powered by wind energy from Iberdrola Australia, the electricity for computations produces negligible carbon waste.
This material is based upon work supported by NSF's LIGO Laboratory, a major facility fully funded by the National Science Foundation. This research has used data, software, and web tools obtained from the Gravitational Wave Open Science Center (https://www.gwopenscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration, and the Virgo Collaboration. LIGO Laboratory and Advanced LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. The Australian Research Council provided additional support for Advanced LIGO. Virgo is funded, through the European Gravitational Observatory (EGO), the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di  (Hunter 2007, v3.2.0), NumPy (Harris et al. 2020, v1.8.1), SciPy (Virtanen et al. 2020, v1.4.1), pandas (pandas development team 2020, v1.0.2), python (Oliphant 2007;Millman & Aivazis 2011, v3.7).