McFACTS Code: AGN-Driven BBH Merger Simulator

• McFACTS Code is a Monte Carlo-based framework that simulates the dynamical evolution of binary black holes within AGN accretion disks.
• It incorporates critical physical processes such as orbital migration, gas damping, and general relativistic inspiral to model merger events.
• The simulations help constrain AGN disk models, black hole population properties, and observational parameters relevant to gravitational-wave detections.

The McFACTS Code (“Monte Carlo For AGN Channel Testing and Simulation”) is a publicly released, modular population-synthesis framework designed to simulate the dynamical evolution and merger history of compact objects—primarily binary black holes (BBHs)—embedded within Active Galactic Nucleus (AGN) accretion disks. It serves as a numerical testbed for quantifying how black holes (BHs) interact with their gaseous and stellar environments, form binaries, and contribute to the gravitational-wave (GW) source population observed by detectors such as LIGO-Virgo-KAGRA (LVK). The code integrates salient physical processes including migration, disk-binary interactions, general relativistic (GR) inspiral, and allows detailed exploration of AGN and nuclear star cluster parameter spaces.

1. Core Features and Population Synthesis

McFACTS initiates galaxy-scale simulations by planting BHs—drawn from user-specified initial mass functions (IMFs) and spin distributions—at assigned orbital parameters (semi-major axis, eccentricity, inclination, periapsis). Supported IMFs include power-law (Pareto-type) distributions, with the capacity to introduce features such as Gaussian bumps to emulate empirical or theorized BH mass spectra. Spin magnitude and orientation assignments support both isotropic and aligned configurations.

AGN disk environments are specified through parameteric models, e.g., the Sirko & Goodman (SG03) or Thompson–Quataert–Murray (TQM05) models, each providing prescriptions for disk surface density $\Sigma(r)$, scale height $h(r)$, viscosity $\alpha$, and optical thickness $\tau$. The code supports both fixed and SMBH-mass-scaled outer/inner disk radii, enabling cosmologically consistent simulations where properties such as the gravitational radius $R_g = GM_{\rm SMBH}/c^2$ play a central role.

The nuclear star cluster (NSC) is represented by its mass and spatial distribution, which regulate the inflow (“capture”) rate of black holes to the disk midplane.

2. Physical Processes: Orbital Evolution, Binary Formation, and Interactions

The time evolution for each embedded BH includes:

• Gas damping: Orbital eccentricity and inclination are damped over timescales approximated by

$$t_{\rm damp} \sim 0.1\,\mathrm{Myr} \left(\frac{q}{10^{-7}}\right)^{-1} \left(\frac{h/r}{0.03}\right)^4 \left(\frac{\Sigma}{10^5\,\mathrm{kg\, m^{-2}}}\right)^{-1} \left(\frac{a}{10^4\,r_g}\right)^{-1/2}$$

where $q = m_{\rm BH}/M_{\rm SMBH}$, $a$ is the orbital radius.

• Migration: Orbital drift is calculated using torque prescriptions that include standard Type I migration and modifications for disk feedback; for example,

$$\frac{\Gamma_{\rm heat}}{\Gamma_{\rm mig}} \propto \frac{1}{\tau \alpha^{3/2}}$$

with $\Gamma_{\rm mig}$ the nominal migration torque.

• Binary formation: Binaries form when neighboring orbits approach within a Hill radius,

$r_H = a_1\left(\frac{M_1}{3M_\bullet}\right)^{1/3}$

where $a_1$ and $M_1$ are the semi-major axis and mass of the more massive BH, $M_\bullet$ is the SMBH mass. Hierarchical mergers are naturally incorporated as earlier-generation remnants encounter additional mergers.

• Dynamical encounters: Encounters among single BHs and binaries produce energy exchanges, leading to binary hardening or, at high velocities, ionization.
• General relativistic inspiral: In the inner disk, close binaries experience orbital decay due to GW emission, using the Peters (1964) formalism for energy loss.

3. Outputs: BBH Merger Properties and Observational Predictions

For every formed binary, McFACTS tracks key observable parameters:

• Mass ratio $q$ and total mass $M_{\rm BBH}$.
• Effective spin $\chi_{\rm eff}$ (aligned with the binary orbital plane), and in-plane spin $\chi_p$.
• GW observables: frequency and characteristic strain, calculated as

$$\nu_{\rm GW} = \frac{\sqrt{GM_{\rm BBH}}}{\pi a_{\rm BBH}^{3/2}}$$

$$h \propto \sqrt{\frac{32}{5}\frac{G^2}{c^4}\frac{M_{\rm BBH}\mu_{\rm BBH}}{D\,a_{\rm BBH}}}$$

where $D$ is the source distance, $\mu_{\rm BBH}$ the reduced mass.

Each simulation realization outputs merger event lists, time evolution histories, generation counts (“1g”, “2g”, etc.), as well as diagnostic plots (e.g., $q$–$\chi_{\rm eff}$ distributions).

To assess observational relevance, simulated mergers are filtered through GW detectability models, employing single-detector SNR thresholds and incorporating cosmological volume and duty cycle corrections. Key relations used include

$$R_{\rm detection} = T_{\rm obs} \sum_k \left[w_{\rm intrinsic,\,k}\,p_{\rm det}(\theta_k)\right]$$

where $T_{\rm obs}$ is observing time, $p_{\rm det}$ the detection probability, and $w_{\rm intrinsic,\,k}$ the intrinsic weight per simulated event.

4. Parameter Studies and Key Findings

Parameter space studies using McFACTS revealed several robust outcomes:

• BBH Mass and Spin Distributions: Flat IMFs and migration pileups predict “echoes” in the mass spectrum—secondary peaks at multiples of the initial pileup mass. Hierarchical mergers progressively increase remnant spins (second generation $\sim$0.7, third $\sim$0.8, fourth $\sim$0.9), matching high-spin GW observations.
• $(q,\chi_{\rm eff})$ Anti-correlation: Simulations consistently produce an anti-correlation in mass ratio versus effective spin, especially pronounced for hierarchical mergers. Fits yield slopes of $d\chi_{\rm eff}/dq\sim-0.25$ for the overall merger sample and up to $-0.84$ when restricted to higher generations—concordant with LVK inferences.
• Disk and NSC Dependence: High-density, short-lived ($0.5$–$2.5$ Myr) disks (e.g., SG03) favor the BBH merger rates and $(q,\chi_{\rm eff})$ trends observed. The dominant contribution originates from SMBHs in $10^7$–$10^{9.4}\,M_\odot$ hosts, where disks are both sufficiently long-lived and supply a high number of embedded BHs.
• Prograde Formation Preference: Simulations indicate that restricting BBH formation to $>90\%$ prograde orientation relative to the disk angular momentum is required to replicate observed anti-correlation structures; even a modest retrograde fraction flattens this relation.

5. Connection to Gravitational-Wave Observations

The McFACTS framework facilitates direct mapping from synthetic merger populations to GW observable distributions. By integrating the likelihood for individual GW events over the simulated population,

$$P(d_j) = \frac{\sum_k w_{{\rm detection},\,k}\,\mathcal{L}(\lambda_k)}{\sum_k w_{{\rm detection},\,k}}$$

in which $\mathcal{L}(\lambda_k)$ is the likelihood for gravitational-wave event parameters $\lambda_k$ (mass, spin, etc.), McFACTS enables inference on the probability that an event such as GW231123 originated in the AGN channel.

Notable findings include the consistency of high-mass, high-spin events such as GW231123 with mergers of third- and fourth-generation black holes in AGN disks, across a range of physically-motivated BH IMFs and lifetimes. For instance, simulations show that producing such events requires AGN disk lifetimes of at least $0.2$–$0.4$ Myr and the participation of multiple hierarchical generations.

6. Impact on Astrophysical Constraints and AGN/NSC Modeling

By matching the rates and distributions of simulated mergers with GW observation catalogs, McFACTS delivers constraints on:

• Disk structural parameters ($\Sigma, h/r, \tau$) and AGN lifetimes ($0.5$–$2.5$ Myr).
• NSC mass and spatial profiles (e.g., preference for steeper profiles to suppress low-mass BH contributions).
• Feedback torques and migration dynamics, including conditions for migration traps and outward disk torques.
• Fraction of retrograde vs. prograde binaries within the nuclear region.

The simulation outcomes feed back into models of AGN-driven supermassive BH growth and feedback in a $\Lambda$CDM galaxy evolution framework, by constraining the efficiency and timescale by which AGN disks process and merge their compact content.

7. Code Infrastructure and Future Directions

McFACTS is implemented in Python for rapid prototyping, transparency, and extensibility. It executes tens of galaxy realizations (per AGN episode) in tens of seconds on modern computing hardware. Its modular architecture allows iterative incorporation of additional processes—such as star–BH interactions, multi-epoch AGN cycles, inclusion of neutron stars, and updates to feedback dynamics.

A plausible implication is that as GW source catalogs grow and mass-spin statistics improve, McFACTS-derived models will enable hierarchical, population-level inference to further isolate AGN contribution among BBH formation channels and to deepen constraints on the astrophysical environments and initial conditions underlying these catastrophic events.