McFACTS: Monte Carlo AGN Channel Simulation

• McFACTS is a population-synthesis framework that models compact-object dynamics in AGN disks, emphasizing binary black-hole formation, migration, and hierarchical mergers.
• It employs a modular Monte Carlo workflow combining microphysical prescriptions with cosmic population synthesis to simulate thousands of black holes through gas torques and dynamical encounters.
• The framework bridges first-principle disk physics and synthetic-universe embedding, enabling direct comparisons with LVK gravitational-wave observations and electromagnetic counterpart predictions.

McFACTS, short for Monte Carlo For AGN Channel Testing and Simulation, is a public, open-source population-synthesis framework for modeling compact-object dynamics in active galactic nucleus (AGN) disks, with emphasis on stellar-mass black-hole (BH) capture, migration, binary black-hole (BBH) formation, hierarchical mergers, and comparison to LIGO–Virgo–KAGRA (LVK) observations. Across the McFACTS series, the framework develops from a one-zone Monte Carlo code for individual SMBH+NSC+disk environments into a synthetic-universe pipeline that computes intrinsic and detection-weighted BBH populations, single-event likelihoods for specific gravitational-wave events such as GW231123, and bolometric electromagnetic (EM) counterpart predictions for AGN-disk embedded mergers (McKernan et al., 2024, Delfavero et al., 2024, Delfavero et al., 19 Aug 2025, McPike et al., 4 Feb 2026).

1. Development and scope

The McFACTS literature is organized as a sequence of progressively broader modeling efforts. The initial code release established the stochastic SMBH+NSC+disk machinery and focused on qualitative signatures of the AGN channel. Subsequent papers added parameter studies of the $(q,\chi_{\rm eff})$ relation, cosmological population synthesis over a synthetic universe, EM counterpart predictions, and event-level inference for GW231123.

Paper	Main scope	Representative result
"McFACTS I: Testing the LVK AGN channel with Monte Carlo For AGN Channel Testing & Simulation (McFACTS)" (McKernan et al., 2024)	One-zone Monte Carlo code for SMBH+NSC+disk realizations	Migration traps or disk boundaries help IMBH growth; McFACTS produces a $(q,\chi_{\rm eff})$ anti-correlation, $\chi_{\rm p}$ tails, and EMRIs
"McFACTS II: Mass Ratio--Effective Spin Relationship of Black Hole Mergers in the AGN Channel" (Cook et al., 2024)	Parameter study of AGN-disk and NSC assumptions	Dense, moderately short-lived SG disks with $\gamma=2$ and prograde fraction $>90\%$ reproduce $d\chi_{\rm eff}/dq=-0.25\pm0.02$
"McFacts III: Compact binary mergers from AGN disks over an entire synthetic universe" (Delfavero et al., 2024)	Synthetic-universe embedding and detection weighting	The majority of observable BBH mergers are expected to originate in galaxies with SMBH mass between $10^{7}M_{\odot}$ and $10^{9.4}M_{\odot}$
"Prospects for the formation of GW231123 from the AGN channel" (Delfavero et al., 19 Aug 2025)	Single-event likelihood analysis for GW231123	GW231123 is consistent with a dynamical BBH merger from the AGN channel; a $(4g,3g)$ origin is most likely for most segregated IMFs and AGN lifetime choices
"McFACTS IV: Electromagnetic Counterparts to AGN Disk Embedded Binary Black Hole Mergers" (McPike et al., 4 Feb 2026)	Bolometric EM counterpart module	In sufficiently dense disks, mergers with chirp mass $\mathcal{M}\gtrsim40M_\odot$ are highly likely to yield observable EM counterparts

The series therefore spans three coupled inferential levels: microphysical prescriptions for migration, hardening, and merger; population-level predictions for mass, spin, redshift, and generation; and event-level or multi-messenger observables. This progression is central to McFACTS’s identity: it is not only a forward simulator of BBH dynamics in AGN disks, but also a detection- and likelihood-aware comparison framework for GW data (Delfavero et al., 2024, Delfavero et al., 19 Aug 2025).

2. Computational architecture and Monte Carlo workflow

In its initial form, McFACTS is a one-zone code in which each galaxy is defined by an SMBH mass, an AGN disk model, an NSC mass and density profile, and an AGN lifetime $(q,\chi_{\rm eff})$0 split into timesteps of duration $(q,\chi_{\rm eff})$1. Each timestep evolves $(q,\chi_{\rm eff})$2–$(q,\chi_{\rm eff})$3 BH through gas torques, dynamical encounters, binary formation and hardening, mergers with GW recoil, and EMRI capture. The code tracks transitions among single eccentric, single circularized, binary, merged remnant, and EMRI categories, and outputs merger catalogs containing masses, spin observables, orbital radius, and time (McKernan et al., 2024).

In the synthetic-universe formulation, McFACTS consists of three major modules: a Universe-and-AGN sampler, an AGN-Disk BBH-Formation Simulator, and a Post-processing & LVK-Detection/Likelihood Module. The universe sampler slices cosmic history into 100 Myr epochs from $(q,\chi_{\rm eff})$4 to $(q,\chi_{\rm eff})$5, samples the galactic stellar-mass function and metallicity evolution, assigns early- or late-type galaxy labels, infers SMBH and NSC masses from empirical scaling laws, and draws an AGN number density from observational luminosity-function constraints. The disk simulator then instantiates a scaled Sirko–Goodman accretion-disk profile with the pAGN disk-modeler, seeds it with a BH population drawn from a chosen initial mass function (IMF), and advances those BH under prescriptions for gas-mediated migration and dynamical encounters until BBH mergers occur, including hierarchical generations. The post-processing module assigns volumetric, intrinsic, and detection weights, computes matched-filter signal-to-noise ratios, and evaluates single-event likelihoods for observed GW events (Delfavero et al., 19 Aug 2025).

The Monte Carlo workflow is explicit. Step 1 slices cosmic history into 100 Myr epochs up to $(q,\chi_{\rm eff})$6. Step 2 samples galaxies, computes the number density of active AGN, instantiates AGN disks, draws BH seeds, and evolves the BH population until BBH mergers occur. Step 3 computes $(q,\chi_{\rm eff})$7, then $(q,\chi_{\rm eff})$8, then the signal-to-noise ratio $(q,\chi_{\rm eff})$9, the detection probability $\chi_{\rm p}$0, and the detection weight $\chi_{\rm p}$1 for each merger. Step 4 convolves the predicted $\chi_{\rm p}$2 distribution with the LVK likelihood $\chi_{\rm p}$3, and Step 5 outputs intrinsic and detection-weighted populations, Bayesian evidence for events, and parameter distributions by merger generation (Delfavero et al., 19 Aug 2025).

A defining design choice is modularity. Alternate disk models, IMFs, or cosmologies can be swapped into the pipeline with minimal changes to the rest of the framework. This makes McFACTS suitable both for controlled parameter sweeps in idealized galaxies and for cosmological forward modeling tied directly to LVK-detectable populations (Delfavero et al., 19 Aug 2025).

3. Physical ingredients and governing prescriptions

McFACTS couples NSC demographics, AGN-disk structure, BH initial conditions, and binary-evolution prescriptions. In the single-galaxy studies, each realization begins with an NSC of total mass $\chi_{\rm p}$4, a broken-power-law radial density profile, and a fixed ratio $\chi_{\rm p}$5. Stellar-mass BH are drawn from a Pareto mass function between 5 and $\chi_{\rm p}$6, with slope varied between $\chi_{\rm p}$7 and $\chi_{\rm p}$8, while initial BH spin magnitudes are sampled from a zero-mean Gaussian with $\chi_{\rm p}$9 and isotropic spin orientations. Initial eccentricities are drawn uniformly from $\gamma=2$0, with $\gamma=2$1 explored from 0 to 0.7 (Cook et al., 2024).

The disk sector uses either the Sirko & Goodman 2003 model or the Thompson–Quataert–Murray 2005 model, implemented via the pAGN library. In McFACTS I, the user may also supply $\gamma=2$2, $\gamma=2$3, and $\gamma=2$4 directly. Typical parameters include $\gamma=2$5, viscosity $\gamma=2$6, disk radial bounds from $\gamma=2$7 to $\gamma=2$8, and an optional migration-trap radius $\gamma=2$9. In the default single-galaxy SG configuration used for EM calculations, the SMBH mass is fixed at $>90\%$0, the AGN accretes at 10% of Eddington, the outer radius is $>90\%$1, the lifetime is $>90\%$2 Myr, and the migration trap occurs at $>90\%$3 (McKernan et al., 2024, McPike et al., 4 Feb 2026).

Binary formation is triggered when two circular, co-orbital or prograde BH approach within a mutual Hill radius. McFACTS I writes this as

$>90\%$4

while McFACTS II gives

$>90\%$5

Once a binary forms, hardening proceeds through gas damping, gas torques, and GW emission. Prograde BH undergo orbital damping with a characteristic timescale

$>90\%$6

and Type I migration is modeled with

$>90\%$7

Retrograde BH can accrete with reversed sign, flip via gas capture, or evolve under analytic dynamical-friction scalings (McKernan et al., 2024).

The IMF sector becomes more elaborate in the GW231123 analysis. McFACTS implements four distinct BH IMFs: OldIMF, with $>90\%$8 on $>90\%$9 plus a Gaussian bump at $d\chi_{\rm eff}/dq=-0.25\pm0.02$0, $d\chi_{\rm eff}/dq=-0.25\pm0.02$1; RomExtended, with $d\chi_{\rm eff}/dq=-0.25\pm0.02$2 on 10–$d\chi_{\rm eff}/dq=-0.25\pm0.02$3; BumpInjection, with $d\chi_{\rm eff}/dq=-0.25\pm0.02$4 plus a Gaussian injection at $d\chi_{\rm eff}/dq=-0.25\pm0.02$5, $d\chi_{\rm eff}/dq=-0.25\pm0.02$6; and CosmicSegregated, a two-step construction in which COSMIC binary-population synthesis is run over metallicity $d\chi_{\rm eff}/dq=-0.25\pm0.02$7 and the resulting BH are weighted by $d\chi_{\rm eff}/dq=-0.25\pm0.02$8 and by $d\chi_{\rm eff}/dq=-0.25\pm0.02$9 to mimic in-spiral segregation, then integrated over cosmic history under Planck 2015 cosmology (Delfavero et al., 19 Aug 2025).

For spin observables, McFACTS uses the standard LVK-aligned definitions

$10^{7}M_{\odot}$0

and

$10^{7}M_{\odot}$1

with $10^{7}M_{\odot}$2. For merger remnants, McFACTS I used Tichy & Marronetti 2008 with a kick of approximately $10^{7}M_{\odot}$3, whereas the GW231123 and EM-counterpart work calls the precession package to compute remnant mass, remnant spin, and GW recoil from the binary masses, spins, and spin-tilt angles (McKernan et al., 2024, Delfavero et al., 19 Aug 2025).

4. Population-level predictions and characteristic signatures

A central McFACTS result is that AGN disks can generate hierarchical merger chains with distinctive mass and spin structure. In the “all-circular, no-dynamics” case, BBH form immediately and merge near the migration trap at $10^{7}M_{\odot}$4 and near the outer boundary at $10^{7}M_{\odot}$5, producing IMBH up to $10^{7}M_{\odot}$6 within $10^{7}M_{\odot}$7 Myr, with mass growth factors $10^{7}M_{\odot}$8–10 from repeated mergers. When retrograde eccentric orbits are allowed, a $10^{7}M_{\odot}$9 Myr delay appears as gas damps eccentricity; when full dynamics is enabled, strong dynamical heating near $10^{9.4}M_{\odot}$0 suppresses trap-driven merging and shifts BBH births and mergers outward to $10^{9.4}M_{\odot}$1 (McKernan et al., 2024).

The mass distribution is strongly IMF-dependent. For a steep IMF with $10^{9.4}M_{\odot}$2, the 1g–1g merger mass function peaks at $10^{9.4}M_{\odot}$3–$10^{9.4}M_{\odot}$4, with a smaller secondary bump at 40–$10^{9.4}M_{\odot}$5, and an overall rate of $10^{9.4}M_{\odot}$6. For a flatter IMF with $10^{9.4}M_{\odot}$7, mergers shift to higher masses, with a dominant peak at $10^{9.4}M_{\odot}$8–$10^{9.4}M_{\odot}$9, a tertiary bump at $(4g,3g)$0–$(4g,3g)$1, a long hierarchical tail to $(4g,3g)$2, and a higher rate of $(4g,3g)$3 (McKernan et al., 2024). In the synthetic-universe calculation, the total-mass distribution peaks at $(4g,3g)$4–$(4g,3g)$5 for first-generation mergers, while hierarchical mergers extend the high-mass tail to $(4g,3g)$6 (Delfavero et al., 2024).

The most extensively quantified diagnostic is the $(4g,3g)$7 relation. In the dense SG “sg_default” run with $(4g,3g)$8 Myr and $(4g,3g)$9, McFACTS II finds distinct islands in the $\mathcal{M}\gtrsim40M_\odot$0 plane: 1g–1g mergers cluster at $\mathcal{M}\gtrsim40M_\odot$1 and $\mathcal{M}\gtrsim40M_\odot$2 with $\mathcal{M}\gtrsim40M_\odot$3; 2g–1g and 2g–2g occupy lower $\mathcal{M}\gtrsim40M_\odot$4 and higher $\mathcal{M}\gtrsim40M_\odot$5, with $\mathcal{M}\gtrsim40M_\odot$6 and $\mathcal{M}\gtrsim40M_\odot$7 for $\mathcal{M}\gtrsim40M_\odot$8; and $\mathcal{M}\gtrsim40M_\odot$9 mergers push further toward $(q,\chi_{\rm eff})$00–0.4 and $(q,\chi_{\rm eff})$01–0.7 for $(q,\chi_{\rm eff})$02. A line through $(q,\chi_{\rm eff})$03 fitted to all events yields

$(q,\chi_{\rm eff})$04

while restricting to hierarchical events steepens the slope to $(q,\chi_{\rm eff})$05 (Cook et al., 2024).

The parameter survey identifies the conditions under which this anti-correlation is reproduced. Flattening the IMF to $(q,\chi_{\rm eff})$06 steepens the full-sample slope to $(q,\chi_{\rm eff})$07 in SG disks. Shorter SG lifetimes of 0.5 Myr flatten the slope to $(q,\chi_{\rm eff})$08, whereas longer lifetimes of 5 Myr steepen it slightly to $(q,\chi_{\rm eff})$09. A retrograde fraction of 0.1 reduces the slope from $(q,\chi_{\rm eff})$10 to $(q,\chi_{\rm eff})$11, and $(q,\chi_{\rm eff})$12 nearly erases the anti-correlation with $(q,\chi_{\rm eff})$13. The preferred scenario is a dense SG-type disk with surface densities $(q,\chi_{\rm eff})$14, lifetime $(q,\chi_{\rm eff})$15–1 Myr, IMF slope $(q,\chi_{\rm eff})$16, negligible retrograde binary formation, initial eccentricities $(q,\chi_{\rm eff})$17, and initial spin width $(q,\chi_{\rm eff})$18 (Cook et al., 2024).

McFACTS also predicts observables beyond $(q,\chi_{\rm eff})$19. Spheroid encounters can tilt embedded BBH orbital planes and generate a tail of $(q,\chi_{\rm eff})$20 up to $(q,\chi_{\rm eff})$21–0.5, whereas turning off spheroid encounters collapses the distribution to $(q,\chi_{\rm eff})$22. Retrograde BH on initially inclined orbits can evolve into EMRIs, with roughly half of the initially retrograde BH becoming EMRIs in the first $(q,\chi_{\rm eff})$23 yr at $(q,\chi_{\rm eff})$24 in the default model. BBH merging at $(q,\chi_{\rm eff})$25 Hz populate the LVK band with strains $(q,\chi_{\rm eff})$26–$(q,\chi_{\rm eff})$27 at $(q,\chi_{\rm eff})$28, while hard binaries in the outer disk can spend thousands of years in the LISA band at $(q,\chi_{\rm eff})$29–$(q,\chi_{\rm eff})$30 Hz (McKernan et al., 2024).

5. Synthetic-universe embedding, detection weights, and event likelihoods

McFACTS III embeds the AGN-channel calculation in a synthetic universe. Galaxies are drawn from the observationally constrained galactic stellar-mass function in 100 Myr bins from $(q,\chi_{\rm eff})$31 to 2, metallicity is assigned via Madau & Dickinson 2014 evolution, galaxy type is classified using a Gaussian-process fit to Peng et al. 2015, SMBH masses are assigned using

$(q,\chi_{\rm eff})$32

and NSC masses are assigned using separate early- and late-type scaling laws capped at $(q,\chi_{\rm eff})$33. AGN number densities are normalized to the infrared bolometric estimate of Lyon et al. 2024, and the AGN duty cycle is set by the assumed lifetime $(q,\chi_{\rm eff})$34 (Delfavero et al., 2024).

The synthetic-universe calculation yields both intrinsic and detection-weighted rate predictions. In the GW231123 formulation, if $(q,\chi_{\rm eff})$35 is the comoving AGN number density and each simulated AGN contributes $(q,\chi_{\rm eff})$36 mergers per disk lifetime $(q,\chi_{\rm eff})$37, then the volumetric weight of merger $(q,\chi_{\rm eff})$38 is

$(q,\chi_{\rm eff})$39

and the intrinsic rate weight is

$(q,\chi_{\rm eff})$40

Using IMRPhenomPv2 and an O3-era PSD, McFACTS computes a matched-filter single-detector $(q,\chi_{\rm eff})$41, a detection probability

$(q,\chi_{\rm eff})$42

and a detection weight

$(q,\chi_{\rm eff})$43

The detection rate is then

$(q,\chi_{\rm eff})$44

summed over samples with $(q,\chi_{\rm eff})$45 (Delfavero et al., 19 Aug 2025).

These calculations imply that AGN disks can account for a substantial subset of the LVK catalog under specific lifetime assumptions. McFACTS III finds that, if hierarchical mergers from AGN disks account for a substantial part of the LVK population, the current models require an AGN lifetime of 0.5 to 2.5 Myr, with $(q,\chi_{\rm eff})$46 and $(q,\chi_{\rm eff})$47 peaking near $(q,\chi_{\rm eff})$48–2.5 Myr. In the favored SG_scaled run, the predicted O3 detection rate is $(q,\chi_{\rm eff})$49–$(q,\chi_{\rm eff})$50 BBH detections yr$(q,\chi_{\rm eff})$51 above $(q,\chi_{\rm eff})$52 if $(q,\chi_{\rm eff})$53 Myr. The majority of observable BBH mergers are expected to originate in galaxies with SMBH mass between $(q,\chi_{\rm eff})$54 and $(q,\chi_{\rm eff})$55 (Delfavero et al., 2024).

Event-level inference is handled by convolving McFACTS’s detection-weighted population with an event likelihood. For GW231123, McFACTS fits the LVK posterior samples obtained with NRSur7dq4 in the four-dimensional space $(q,\chi_{\rm eff})$56 using a truncated multivariate Gaussian, and computes

$(q,\chi_{\rm eff})$57

in practice as

$(q,\chi_{\rm eff})$58

Bayes factors between two McFACTS models are then formed as $(q,\chi_{\rm eff})$59. Using this machinery, the GW231123 study finds that the event is consistent with a dynamical BBH merger from the AGN channel and postulates that its masses and spin magnitudes are most consistent with a merger of fourth- and third-generation BH for most choices of a segregated BH IMF and AGN lifetime (Delfavero et al., 19 Aug 2025).

McFACTS III also places earlier high-mass LVK events in this framework. Many O3 events lie within the 68% contours of the predicted population in the $(q,\chi_{\rm eff})$60 plane, including GW190929_012149, which is naturally explained by a 1g–2g or 2g–2g merger in an AGN disk, and GW190521, which lies at the edge of the intrinsic high-mass tail but is still produced at low probability by hierarchical mergers (Delfavero et al., 2024).

6. Electromagnetic extension, assumptions, and limitations

McFACTS IV extends the framework to predict bolometric EM luminosities from AGN-disk embedded BBH mergers. Two channels are modeled. The first is recoil-driven shock luminosity: a GW recoil kick perturbs the merger remnant’s orbit through the disk and drives a ram-pressure shock. The energy and timescale are estimated as

$(q,\chi_{\rm eff})$61

and

$(q,\chi_{\rm eff})$62

giving

$(q,\chi_{\rm eff})$63

The second is jet luminosity from high-spin remnants accreting from the disk. McFACTS estimates a Bondi–Hoyle–Lyttleton rate

$(q,\chi_{\rm eff})$64

and adopts

$(q,\chi_{\rm eff})$65

with $(q,\chi_{\rm eff})$66, $(q,\chi_{\rm eff})$67–3 for $(q,\chi_{\rm eff})$68, and a conservative jet lifetime of $(q,\chi_{\rm eff})$69 s (McPike et al., 4 Feb 2026).

Jet observability depends on breakout through the disk. McFACTS uses

$(q,\chi_{\rm eff})$70

with $(q,\chi_{\rm eff})$71, and a breakout time $(q,\chi_{\rm eff})$72. If $(q,\chi_{\rm eff})$73, the jet can be observed; otherwise the energy diffuses out on $(q,\chi_{\rm eff})$74, with a diffusion-limited luminosity $(q,\chi_{\rm eff})$75 (McPike et al., 4 Feb 2026).

The resulting EM statistics are highly mass dependent. Shock luminosities are typically $(q,\chi_{\rm eff})$76, negligible compared to AGN variability, while jet luminosities span $(q,\chi_{\rm eff})$77–$(q,\chi_{\rm eff})$78, with the highest-generation remnants near the migration trap producing $(q,\chi_{\rm eff})$79. Under the default SG model, the probability of an observable EM flare rises from $(q,\chi_{\rm eff})$80–$(q,\chi_{\rm eff})$81 to $(q,\chi_{\rm eff})$82–$(q,\chi_{\rm eff})$83, $(q,\chi_{\rm eff})$84–$(q,\chi_{\rm eff})$85, and $(q,\chi_{\rm eff})$86–$(q,\chi_{\rm eff})$87. McFACTS IV therefore states that BBH with $(q,\chi_{\rm eff})$88 are essentially guaranteed to power observable jets, neglecting viewing-angle and Type-II AGN obscuration, and recommends daily–weekly monitoring for up to $(q,\chi_{\rm eff})$89 days post-merger in surveys such as ZTF or LSST (McPike et al., 4 Feb 2026).

The framework’s limitations are explicit. In the GW231123 likelihood analysis, all AGN disks are assumed to follow the same scaled Sirko–Goodman profile and share a single lifetime $(q,\chi_{\rm eff})$90; BH spin-tilt angles are not tracked, and spin magnitudes alone are used because GW231123’s posteriors for $(q,\chi_{\rm eff})$91 versus $(q,\chi_{\rm eff})$92 disagree among waveform models; and the approximate nature of $(q,\chi_{\rm eff})$93 and of frozen-mass $(q,\chi_{\rm eff})$94 limits precise absolute rates, though relative-rate trends are robust. The hierarchical merger prescriptions assume isotropic spin orientations upon binary formation inside the disk, aside from preserving the in-plane coherence imparted by gas torques. No explicit modeling of disk self-gravity or of feedback from repeated mergers on disk structure is included, and the CosmicSegregated IMF inherits uncertainties from COSMIC’s binary-evolution assumptions. In the EM module, jets and shocks are treated in a one-zone, bolometric manner without full radiative-transfer or spectral modeling, and AGN variability, viewing-angle effects, Type-II obscuration, and relativistic beaming are not rigorously included (Delfavero et al., 19 Aug 2025, McPike et al., 4 Feb 2026).

Taken together, McFACTS provides a self-consistent bridge from first-principle IMF and disk-physics assumptions, through Monte Carlo–driven hierarchical merger simulations, to LVK-detectable BBH populations, single-event likelihoods, and, in its latest extension, joint GW–EM phenomenology. This suggests its primary scientific role is not merely catalog generation, but parameter inference on BH mass segregation, AGN episodic lifetimes, disk structure, and the role of gas in dynamical BBH assembly (Delfavero et al., 19 Aug 2025).