The paper
BAQARO: Tracing Stochastic Black Hole Growth Histories and Quasar Lightcurves in a Cosmological Context
Elia Pizzati, Joseph F. Hennawi & Joop Schaye. MNRAS, submitted.
Black hole Accretion and Quasar Activity in a Realistic Observational framework: a semi-empirical model for how supermassive black holes grow and light up as quasars, from cosmic dawn to the present day.
BAQARO is built on the subhalo merger trees of the FLAMINGO-10k dark-matter-only simulation, a trillion particles in a 2.8 Gpc box, large enough to contain the rare, massive haloes that host the brightest quasars. A black hole is seeded in each halo as it enters the tree, and grows through accretion tied to its host’s cold-gas supply. Six free parameters are fit jointly to the bolometric quasar luminosity function, quasar clustering, and the conditional Eddington-ratio distribution at fixed luminosity, over 0 ≲ z ≲ 7.
The defining choice is that accretion is treated as intrinsically stochastic: rather than smooth, deterministic growth tracks, black holes accrete in fluctuating episodes set by a coherence timescale and an overall scatter. That is not a detail. It is what lets rare billion-solar-mass black holes emerge as outliers of the accretion-rate distribution while the bulk of the population grows far more steadily. Those same fluctuations set the quasar lightcurve, so one timescale governs both how a black hole grows and how long it shines, leaving its imprint on lightcurves, on the proximity zones a quasar carves into the intergalactic medium, and on quasar clustering.
Mergers are followed alongside accretion, on the simulation's own trees. They stay subdominant, mattering most at z ≲ 1 for the heaviest black holes once cold gas has run out. Tracking them also produces a catalogue of coalescing black hole binaries, connecting the model directly to the gravitational-wave channel.
The code
The full pipeline: forward model, Gaussian-process emulators, likelihoods and MCMC, and the scripts that make every figure in the paper.
Browse on GitHub →How the model works
Halos in, quasars out
What goes in is a subhalo merger tree, and nothing else: no galaxies, no gas
physics, no tuning against redshift. What comes out is a black hole in every
subhalo, carried forward on timesteps set by the coherence time rather than by
the simulation's snapshots, with a mass, an accretion rate and a luminosity
at each one, and a lightcurve joining them up. The figure above is the whole thing at a
glance: on the left, the FLAMINGO-10k merger trees, the seeding, and
the stochastic accretion that sets both the mass growth and the lightcurve;
top right, the six parameters, emulated and compared against the
luminosity function, quasar clustering and the conditional ERDF;
bottom right, what the fitted model then predicts.
Between subhalo and quasar sit three ingredients, and six numbers.
Seeding sets the mass a black hole is born with, through a seed mass
Mseed and its scatter σseed.
Accretion fixes how fast it then grows. The accretion rate
ηacc is drawn from a lognormal whose mean is a power law
in the host's specific cold accretion rate sṀcold:
ηav,0 is that mean at a reference supply of
1 Gyr−1, ηav,slope is the
power-law index, and σacc is the scatter about it.
The coherence time τcoherence fixes how long a single accretion episode
lasts before the rate is redrawn. Those six are the entire model, and they
are what the fit determines.
The sliders below are there to help visualise how each parameter is defined
and what it changes.
1 Seeding: the mass a black hole is born with
Every subhalo gets a black hole when it first appears in the tree, with a
mass a fixed fraction of the subhalo's own and lognormal scatter about it.
Subhalos are identified close to the resolution threshold, so the width of
the curve is essentially that scatter: log Mseed
slides it left and right, quoted at a pivot subhalo mass of
1010.5 M☉, and
σseed sets its width.
This is an empirical prescription. The simulation does not resolve
the scales on which seeds actually form, so Mseed is not a
physical seed mass but an anchor marking where a growth track starts. The
second curve is the mass function at z = 0, scaled to the
same total number density: the same black holes, 13 Gyr later.
2 Accretion: how fast it grows
A black hole accretes at a rate ηacc = ṀBH,acc/ṀEdd,
drawn from a lognormal whose mean follows the specific cold accretion rate
sṀcold of its host, so black holes grow fastest when their
hosts do. The cold fraction
fcold(Mhalo, z) is taken
from Correa et al.
(2018): in massive haloes the inflow shock-heats to the virial
temperature and settles into a hot atmosphere, so less of it ever reaches
the centre cold.
This relation between ηacc and sṀcold
is the same at every redshift. What evolves is the supply:
sṀcold falls as the universe ages, and the points in the
plot below trace its population median from z = 8 down to
z = 0.
log ηav,0 raises or lowers the relation,
ηav,slope tilts it, and
σacc is its width.
3 Coherence time: how long each accretion episode lasts
The accretion rate is redrawn on the coherence time τ, not once
per simulation snapshot. Snapshots are tens to hundreds of Myr apart and
τ is far shorter, so a black hole gets many independent draws
between them. The result resembles a real quasar lightcurve: a sequence of
coherent bursts. In reality many timescales contribute at once, and the
approximation here is to let a single one drive the process. Each flat segment is one
draw, its length set by τ, its level by
ηav,0 and ηav,slope, and its
scatter by σacc.
Accretion rate and luminosity are not the same thing. ηacc is the rate; λEdd is the light it produces. Radiative efficiency depends on the rate, so the two separate most in the brightest bursts, which is how the model grows black holes fast without over-producing luminous quasars.
The third curve is ηacc averaged from the start up to each moment; the dashed line is what that average tends to over a long enough stretch. τ controls how quickly one reaches the other. Short τ packs hundreds of draws into the window, so the average settles at once and every black hole grows alike. Long τ leaves a handful, so it can sit well above the line for all 100 Myr. That spread is what lets an individual black hole end up far above the average, and it is why τ shapes the massive end of the population.
The lower panel follows the consequences, for a black hole starting at 108 M☉. Its mass integrates the accretion rate, so time above the average is banked and never lost: the one that opened with two strong bursts stays ahead. Beside it is the lightcurve the same history produces, which follows the rate directly instead of integrating it, so the bursts that build the mass are also when the object shines as a bright quasar.
Mergers are tracked along the trees as well, but they turn out to contribute only marginally to black hole growth compared to accretion, even under generous assumptions about merger timescales and remnant survival.
Observable predictions
The three panels above show the ingredients of the model. The panels below show what those ingredients produce when the model is run over the whole halo catalogue: the summary statistics that the fit is made against. Dashed lines mark the best-fit model throughout.
Running the model is expensive: the fiducial run writes 780 GB across sixteen nodes, and even a subsampled run takes hours. The summary statistics are therefore emulated with Gaussian processes trained on a few thousand forward runs, which reproduce the model in milliseconds. This is what makes the inference below possible. It is also what powers these sliders: move any of the six parameters and every panel is recomputed in your browser. Held out from the training set, the emulators reproduce the forward model to 0.11 dex (luminosity function), 0.22 dex (Eddington ratios) and 0.20 dex (host halos) root-mean-square — smaller than the observational uncertainties they are fitted against.
What it predicts
Quasar luminosity function
redshift (up to 7)
threshold
L > 10^45.5 erg/s
Fitting the model
The model is fitted to three datasets at once: the quasar luminosity function (QLF), the conditional Eddington-ratio distribution (cERDF) and quasar clustering (Corr), the same three quantities shown in the sliders above. The posterior is sampled for every combination of them, so the constraining power of each combination of observables can be assessed.
The fiducial model
The fiducial model is the best-fit combination from the joint fit to all three (yellow stars in the corner plot). It is what the sliders return to, what the panels above are drawn at, and what every result below comes from.
| Parameter | Best fit | Prior | Meaning |
|---|---|---|---|
| log10(Mseed/M☉) | 4.03 | [3.0, 7.0] | seed mass at the pivot halo mass 10^10.5 M_sun |
| σseed | 0.51 | [0.2, 1.5] | scatter in the seed mass at fixed halo mass [dex] |
| log10 ηav,0 | -1.24 | [-2.5, -0.5] | average accretion rate at a cold-gas supply of 1/Gyr |
| ηav,slope | 0.83 | [0.2, 2.0] | how that average tracks the cold-gas supply |
| σacc | 0.51 | [0.2, 1.0] | scatter in the accretion rate [dex] |
| log10(τcoherence/yr) | 5.90 | [3.5, 7.3] | coherence timescale of accretion |
Results
Fitted to the luminosity function, the Eddington-ratio distribution and clustering, the model reproduces all three from z ≈ 7 to z ≈ 0. It also reproduces quantities it was never fitted to: the local black hole mass function, the local scaling relations, and the cosmic black hole mass and accretion densities. Three conclusions stand out.
- Growth follows the halos. Accretion tracks the cold gas arriving in halos, with no explicit dependence on cosmic time, so the redshift evolution of quasars, downsizing included, emerges from the halos and the cold gas they supply.
- Massive black holes are early winners of a lottery. A 109 M☉ black hole at z ≈ 6 comes from a lucky sequence of short, radiatively inefficient super-Eddington bursts lasting of order 1 Myr, not from sustained Eddington-limited growth. The same coherence time shapes the quasar lightcurves, and leaves its imprint on proximity zones and clustering.
- Mergers are a minor channel. They contribute little to growth except at z ≲ 1 for the most massive black holes, once cold gas has run out. They do, however, contribute to shaping the high-mass tail of the black hole mass function.
Not everything matches. The model over-produces the faint end of the luminosity function relative to pre-JWST surveys, although JWST may be finding exactly the faint AGN those surveys missed. Its Eddington ratios sit slightly too high at z ≈ 5 (where the measured values rest on C IV virial masses) and too low at z ≈ 1, where their distribution is also broader than observed, especially for bright sources. And the clustering of z ≈ 4 quasars remains somewhat stronger than the model predicts. Each of these may reflect observational systematics as much as a limitation of the model, and will need further work on both the observational and the modelling side.
Data products coming soon
Summary statistics on the model's own grids, the portable emulators behind the sliders above, the MCMC chains, and a per-object catalogue of the bright quasars across 41 redshifts. Every file carries a provenance record of the parameters, code revision and build time that produced it. These are released once the paper is accepted, and available on request before then.
Citation
If you use BAQARO, please cite the model paper. If you use the data products, please cite its DOI as well.