pith. sign in

arxiv: 2605.24112 · v1 · pith:CRC4ZUEWnew · submitted 2026-05-22 · 🌌 astro-ph.GA · astro-ph.CO· astro-ph.HE

The Lumina Project: The Demographics of Active Galactic Nuclei from Quasars to Little Red Dots at zgeq 3

Xuejian Shen , Oliver Zier , Aaron Smith , Rongrong Liu , Rahul Kannan , Teodora-Elena Bulichi , Sonja M. Koehler , Volker Springel
show 11 more authors
Mark Vogelsberger Lars Hernquist Rohan P. Naidu Anna de Graaff Elia Pizzati David M. Alexander Luis C. Ho Vasily Kokorev Gene Leung Anna-Christina Eilers Ryan C. Hickox
This is my paper

Pith reviewed 2026-06-30 15:10 UTC · model grok-4.3

classification 🌌 astro-ph.GA astro-ph.COastro-ph.HE
keywords active galactic nucleilittle red dotsluminosity functionshigh-redshiftsupermassive black holescosmological simulationsreionization
0
0 comments X

The pith

An empirical model from a large cosmological simulation reproduces AGN luminosity functions for both quasars and Little Red Dots at z ≥ 3.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper builds an empirical mapping from simulated supermassive black holes to observed active galactic nuclei populations. Low-mass black holes are assigned to a Little Red Dot phase with a 30 percent duty cycle, combined with bolometric corrections and a 0.3 dex scatter in luminosity. This framework matches the observed luminosity functions and clustering of Little Red Dots while recovering earlier quasar constraints. It also tracks the contribution of these objects to black hole mass density and helium reionization. A reader would care because it offers a single population that connects faint JWST sources to brighter pre-JWST quasars within one cosmological volume.

Core claim

The central claim is that SMBHs with masses up to 10^7 solar masses remain in the LRD phase with a 30 percent duty cycle, and that applying standard bolometric and extinction corrections plus 0.3 dex log-normal scatter in bolometric luminosity allows the simulated population to reproduce the observed AGN luminosity functions and clustering from quasars down to Little Red Dots at z ≥ 3, while the same population drives the modeled He II reionization.

What carries the argument

The empirical AGN model that assigns SMBHs below 10^7 solar masses to an LRD phase with 30 percent duty cycle and applies luminosity corrections plus scatter to map simulated black holes onto multi-band observations.

If this is right

  • The modeled AGN population is the dominant driver of He II reionization in the simulation.
  • Number densities of AGN and LRDs evolve with redshift in a way that can be directly compared to future surveys.
  • The integrated black hole mass density receives a calculable contribution from the LRD phase at high redshift.
  • The same framework supports building general population-synthesis models of high-redshift AGN that include LRDs.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The model implies that Little Red Dots mark a distinct early growth stage for black holes that later transition into quasars.
  • Varying the duty cycle with redshift or mass could alter the timing of helium reionization in testable ways.
  • The required luminosity scatter suggests that future multi-epoch observations might detect variability differences between LRDs and quasars.
  • Extending the model to include feedback effects on host galaxies could link AGN demographics to galaxy quenching at these epochs.

Load-bearing premise

Supermassive black holes with masses up to 10 million solar masses remain in the Little Red Dot phase with a fixed 30 percent duty cycle.

What would settle it

An observation that the spatial clustering amplitude of Little Red Dots at z ≥ 3 deviates significantly from the prediction based on the underlying galaxy distribution in the simulation, or that the bright-end quasar luminosity function cannot be recovered with 0.3 dex scatter.

Figures

Figures reproduced from arXiv: 2605.24112 by Aaron Smith, Anna-Christina Eilers, Anna de Graaff, David M. Alexander, Elia Pizzati, Gene Leung, Lars Hernquist, Luis C. Ho, Mark Vogelsberger, Oliver Zier, Rahul Kannan, Rohan P. Naidu, Rongrong Liu, Ryan C. Hickox, Sonja M. Koehler, Teodora-Elena Bulichi, Vasily Kokorev, Volker Springel, Xuejian Shen.

Figure 1
Figure 1. Figure 1: DM surface density field of a slice of the simulation volume centred on the most luminous AGN at 𝑧 ≃ 4. The thickness of the slice is 150 cMpc. We first zoom into a 75 × 75 cMpc2 subfield and present the blending of the mass-weighted He iii fraction and the gas surface density field. This illustrates the He iii bubble driven by the hard ionizing radiation from the central AGN. We further zoom in on the 25 … view at source ↗
Figure 2
Figure 2. Figure 2: A blending of the gas surface density map and the mass-weighted fraction of He iii for a slice of the simulation volume centred on the most luminous AGN at 𝑧 ≃ 4. The He iii fraction is colour-mapped while the gas surface density dictates the transparency. The thickness of the slice is 150 cMpc. This image illustrates the large-scale topology of He ii reionization driven mostly by AGN in Lumina. et al. 201… view at source ↗
Figure 3
Figure 3. Figure 3: BH lightcone wedges from Lumina. The wedge shows a 3.6 ◦ × 3.6 ◦ pencil-beam (full opening angle 0.063 rad) cut from the lightcone (corrected for redshift space distortions) from 𝑧 ≃ 3 to 𝑧 ≃ 6. Each filled circle marks a BH with sightline comoving distance plotted along the long axis and the transverse position along the short axis. Marker colour encodes the bolometric luminosity, while marker size encode… view at source ↗
Figure 4
Figure 4. Figure 4 [PITH_FULL_IMAGE:figures/full_fig_p009_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Left: Hydrogen column density distribution function at 𝑧 ≃ 4 at intrinsic hard X-ray luminosity 𝐿 int HX = 1043.5 erg s−1 . The distribution is normalised over log (𝑁H [cm−2 ] ) = [20, 26]. We show the 𝑁H distribution from Ueda et al. (2003) (used in Hopkins et al. 2007), Ueda et al. (2014) (used in Shen et al. 2020a), and more recent constraints from Ananna et al. (2019) leveraging the NuSTAR data. We sho… view at source ↗
Figure 6
Figure 6. Figure 6: A schematic view of the empirical model for AGN LF. We start from the raw intrinsic bolometric luminosities of SMBHs informed from Lumina and apply a scatter in bolometric luminosities due to unresolved variability of accretion (see Section 3.1). Then we perform a split of canonical AGN and the ones in the LRD phase based on the BH particle mass in Lumina. The former undergoes the standard AGN bolometric a… view at source ↗
Figure 7
Figure 7. Figure 7: SMBH mass function at 𝑧 ≃ 5 in Lumina, TNG100, and Mtng. We compare the results with observational constraints of BLAGN from Wu et al. (2022) for quasars and Matthee et al. (2024); Taylor et al. (2025) for LRDs. The latter should be considered as upper limits due to the large uncertainties in single-epoch BH mass estimates using local scaling relations. Specifically, we highlight the upper mass threshold (… view at source ↗
Figure 8
Figure 8. Figure 8: Bolometric LFs of AGN at 𝑧 = 3 − 6. We present the bolometric LFs of SMBHs in Lumina, TNG100, and Mtng. The blue dashed line shows the contribution by LRDs. These simulation predictions are compared to the best-fit quasar LF in Shen et al. (2020a) (the “Global fit A”) and the compiled observational data there (moved to the bolometric plane by preserving the offsets to the best-fit model). We also compare o… view at source ↗
Figure 9
Figure 9. Figure 9: Rest-frame optical (B band, ∼ 4400Å) LFs of AGN at 𝑧 ≃ 3 − 5 from Lumina, TNG100, and Mtng. The shaded region represents the lower and upper limits of the LFs given uncertainties in bolometric corrections of the LRDs. We compare the results to optical observations compiled in Shen et al. (2020a) and LRD constraints from Kokorev et al. (2024); Ma et al. (2025) (at 5100Å) and Matthee et al. (2024); Lin et al… view at source ↗
Figure 10
Figure 10. Figure 10: Rest-frame UV LFs of AGN at 𝑧 = 3 − 6 from Lumina, TNG100, and Mtng. We compare them with the UV observations compiled in Shen et al. (2020a) and additional constraints from Niida et al. (2020); Schindler et al. (2023). For reference, we also show the best-fit AGN UV LF from Shen et al. (2020a) and Kulkarni et al. (2019). For LRDs, we show two groups of observational constraints: the BLAGN LFs from Matthe… view at source ↗
Figure 11
Figure 11. Figure 11: The luminosity ratio between BH and host galaxy in the UV versus in the optical B band at 𝑧 ≃ 5. We show the distribution of BHs below the LRD selection mass threshold 𝑀crit BH on this plane. The majority of the BHs are dominated by their host galaxies in the UV and optical, which raises the concern that the fraction of sources that would pass observational selection criteria of LRDs is smaller than the 𝑓… view at source ↗
Figure 12
Figure 12. Figure 12: Hard X-ray (2 − 10 keV) LFs of AGN at 𝑧 = 3 − 6 from Lumina, TNG100, and Mtng. We note that the X-ray LFs here are only contributed by the non-LRD population we identify in the simulations, so there is no dependence on the 𝑓duty and bolometric corrections of LRDs. We compare the simulation results with X-ray observations compiled in Shen et al. (2020a) along with more recent constraints from Ananna et al.… view at source ↗
Figure 13
Figure 13. Figure 13: Rest-frame mid-IR band (∼ 15𝜇m) LFs of AGN at 𝑧 = 3 − 5 in Lumina. Results from TNG100 and Mtng are shown for reference. For com￾parison, we include observational constraints compiled in Shen et al. (2020a). We also include JWST MIRI constraints from Bulichi et al. (2026), with data converted from their bolometric LFs or from their 3𝜇m LFs assuming the canonical AGN SED (Krawczyk et al. 2013). 4.5 Mid-IR … view at source ↗
Figure 14
Figure 14. Figure 14: Projected cross-correlation function (𝑤p) of LRDs with galaxies at 𝑧 ≃ 4 (left panel) and 𝑧 ≃ 5 (right panel) in Lumina. As discussed in the main text, we select galaxies with SFR> 1.7 M⊙ yr−1 to mimic the tracer galaxy sample used in observations. We further measure the bias factor of the selected galaxy sample and apply a 𝑏 sim gal /𝑏 obs gal correction to the observed 𝑤p in Lin et al. (2026b) and Arita… view at source ↗
Figure 15
Figure 15. Figure 15: Number density of AGN versus redshift. We show the number den￾sity of all AGN with 𝐿bol ≥ 1043 erg s−1 , and LRDs with different thresholds of 𝐿bol as labelled. For comparison, we show observational constraints for LRDs from Kokorev et al. (2024); Greene et al. (2026) at 𝑧 ≳ 5 and Ma et al. (2026); Zhuang et al. (2026) at lower redshifts. The constraints from Kokorev et al. (2024); Greene et al. (2026) re… view at source ↗
Figure 16
Figure 16. Figure 16: SMBH accretion rate density (BHAD) versus redshift, obtained through a summation of accretion rates of BHs in Lumina in the bolometric luminosity range 𝐿bol = 1043 −1048 erg s−1 . We show the BHAD contributed by all AGN and by LRDs in Lumina and compare them to LRD constraints compiled in Inayoshi & Ichikawa (2024) (originally taken from Matthee et al. 2024; Greene et al. 2024; Kokorev et al. 2024; Akins … view at source ↗
Figure 18
Figure 18. Figure 18: Comoving production rate density of He ii-ionizing photons at 54.42 − 100 eV, 𝑛¤ HeII ion , as a function of redshift in Lumina. We show decom￾posed emissivity from SMBHs in four mass bins. The vertical shaded band marks the He ii-reionization epoch in Lumina, 𝑧 ≃ 3 − 6. We also show the contribution from other sources to 𝑧 ≃ 5, which are subdominant for He ii reionization. They are therefore not tracked … view at source ↗
Figure 19
Figure 19. Figure 19: Bolometric LFs of AGN at 𝑧 ≃ 3 in Lumina, assuming 𝜎bol = 0 and 0.7 dex, compared to the fiducial model. We choose 𝑧 ≃ 3 for this comparison because it is the redshift where Lumina predictions have most overlap with the bright-end observational constraints. A lower (higher) 𝜎bol results in under- (over-) prediction of the bright-end LF, while the faint end is less affected. A luminosity- or mass-dependent… view at source ↗
Figure 20
Figure 20. Figure 20: Top: AGN B-band optical LF at 𝑧 ≃ 5 in Lumina, for several combinations of 𝑓duty and 𝜅. For clarity, we show only the predictions using the “weak IR” bolometric correction of LRDs. The observational constraints of LRDs compiled in Section 4.2 are shown here for comparison. 𝑓duty controls the normalization of the faint-end LF that is dominated by LRDs, while 𝜅 mainly influences the transitional luminosity … view at source ↗
Figure 21
Figure 21. Figure 21: Comparison of AGN bolometric LFs at 𝑧 ≃ 4 and 6 between Lumina and ASTRID (Ni et al. 2022, 2025), FLAMINGO (Schaye et al. 2023; Ding et al. 2026), BRAHMA (Bhowmick et al. 2024b), and FIRE-2 (a “passive” BH model; Marszewski et al. 2026). We include the observational data compiled in Section 4.1 for comparison. combined with a self-consistent treatment of super-Eddington phases (as we discussed in Section … view at source ↗
read the original abstract

High-redshift active galactic nuclei (AGN) serve as powerful probes of early black-hole growth, galaxy formation, and the evolving intergalactic medium (IGM). In this work, we use Lumina, a cosmological radiation-hydrodynamic simulation spanning the epochs of hydrogen and helium reionization, which combines a large $(500\,{\rm cMpc})^3$ volume with $2\times 6000^3$ resolution elements, to explore high-redshift AGN. The simulation self-consistently follows hundreds of millions of galaxies and supermassive black holes (SMBHs), together with their impact on the ionization and thermal state of the IGM. We exploit this uniquely large dynamic range to predict multi-band AGN luminosity functions (LFs) at $z \geq 3$, from hard X-rays to the mid-infrared. These predictions encompass both moderately luminous quasars and the faint ``Little Red Dots'' (LRDs) uncovered by JWST. We develop an empirical model that maps simulated SMBHs onto observed AGN using bolometric and extinction/absorption corrections for canonical AGN and LRDs, and in which SMBHs with $M_{\rm BH}\leq 10\,M_{\rm seed} \sim 10^{7}\,{\rm M}_{\odot}$ stay in the LRD phase with a duty cycle of $30\%$. This simple framework reproduces the observed LFs and clustering of LRDs. Meanwhile, the pre-JWST quasar LF constraints are recovered, although we find that a $\sim 0.3$ dex log-normal scatter in bolometric luminosity is required to reproduce the bright end. We place the simulated AGN population in the cosmological context by quantifying the redshift evolution of AGN and LRD number densities, and their contributions to the integrated BH mass densities. The same AGN population is the dominant driver for the HeII reionization modelled self-consistently in Lumina. This empirical AGN model paves the way for general population-synthesis models of high-redshift AGN, including LRDs, in a unified cosmological framework.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

3 major / 1 minor

Summary. The paper presents results from the Lumina cosmological radiation-hydrodynamic simulation (500 cMpc volume, 2×6000³ resolution) to predict multi-band AGN luminosity functions at z≥3, encompassing both quasars and JWST-discovered Little Red Dots (LRDs). An empirical model maps simulated SMBHs to observed AGN via bolometric/extinction corrections, with SMBHs below ~10^7 M_⊙ assigned a 30% duty cycle in the LRD phase; a ~0.3 dex log-normal scatter in bolometric luminosity is introduced to match the bright end. The framework is stated to reproduce observed LFs and LRD clustering, recover pre-JWST quasar constraints, quantify AGN number densities and BH mass density contributions, and identify the simulated AGN as the dominant driver of self-consistently modeled HeII reionization.

Significance. If the empirical mapping can be shown to have predictive power beyond parameter tuning, the work supplies a unified cosmological framework linking high-z quasars and LRDs within a single large-volume simulation that self-consistently treats reionization. The large dynamic range and explicit connection to IGM thermal/ionization state are clear strengths for population-synthesis modeling.

major comments (3)
  1. [Abstract / empirical model] Abstract and empirical model description: the 30% duty cycle for SMBHs with M_BH ≤ 10 M_seed ~10^7 M_⊙ is stated to be chosen so that the model reproduces the JWST LRD luminosity function. Because the central claim is that the framework 'reproduces the observed LFs', this choice renders the match partly by construction; an independent validation (e.g., using the duty cycle fixed from other observables and then predicting the LF) is required for the reproduction to be non-circular.
  2. [Abstract] Abstract: the statement that '~0.3 dex log-normal scatter in bolometric luminosity is required to reproduce the bright end' is load-bearing for the quasar LF recovery claim, yet no quantitative comparison (with vs. without scatter) or goodness-of-fit metric is referenced. This needs to be shown explicitly, for example via a table or figure comparing model LFs at the bright end.
  3. [Abstract] Abstract: the claim that clustering of LRDs is reproduced must be shown to be independent of the duty-cycle tuning. If the duty cycle only rescales number density without altering the spatial bias of the selected objects, this should be demonstrated; otherwise any additional parameter adjustment for clustering statistics must be reported.
minor comments (1)
  1. Notation for M_seed and the precise definition of the LRD phase should be clarified in the main text to avoid ambiguity with the abstract's 'M_BH ≤ 10 M_seed ~10^7 M_⊙'.

Simulated Author's Rebuttal

3 responses · 0 unresolved

We thank the referee for their careful reading and valuable comments on our manuscript. We address each of the major comments below and indicate the revisions we will make.

read point-by-point responses

  1. Referee: [Abstract / empirical model] Abstract and empirical model description: the 30% duty cycle for SMBHs with M_BH ≤ 10 M_seed ~10^7 M_⊙ is stated to be chosen so that the model reproduces the JWST LRD luminosity function. Because the central claim is that the framework 'reproduces the observed LFs', this choice renders the match partly by construction; an independent validation (e.g., using the duty cycle fixed from other observables and then predicting the LF) is required for the reproduction to be non-circular.

    Authors: We agree that the 30% duty cycle parameter was chosen specifically to reproduce the observed JWST LRD luminosity function, as explicitly stated in the abstract and methods. This does make the match to the LRD LF partly by construction. The framework's strength lies in its ability to simultaneously match multiple observables, including the quasar LF, clustering, and reionization history, within a single self-consistent simulation. For the revision, we will add a new subsection discussing independent motivations for the duty cycle value (e.g., from AGN variability studies at lower redshifts) and present results for a range of duty cycle values to show robustness. revision: yes

  2. Referee: [Abstract] Abstract: the statement that '~0.3 dex log-normal scatter in bolometric luminosity is required to reproduce the bright end' is load-bearing for the quasar LF recovery claim, yet no quantitative comparison (with vs. without scatter) or goodness-of-fit metric is referenced. This needs to be shown explicitly, for example via a table or figure comparing model LFs at the bright end.

    Authors: We acknowledge that the manuscript does not include an explicit quantitative comparison of the luminosity function with and without the scatter. In the revised manuscript, we will add a figure showing the bright-end quasar LF both with and without the 0.3 dex log-normal scatter, along with a table of goodness-of-fit metrics (e.g., reduced chi-squared) to the observational data points. revision: yes

  3. Referee: [Abstract] Abstract: the claim that clustering of LRDs is reproduced must be shown to be independent of the duty-cycle tuning. If the duty cycle only rescales number density without altering the spatial bias of the selected objects, this should be demonstrated; otherwise any additional parameter adjustment for clustering statistics must be reported.

    Authors: The duty cycle is a uniform rescaling factor applied to the number of active SMBHs and does not alter the spatial distribution or halo bias of the underlying population, which is set by the simulation. We will demonstrate this explicitly in the revised version by computing the two-point correlation function for the LRD population at different duty cycle values (e.g., 10%, 30%, 50%) and showing that the bias factor remains consistent within uncertainties. revision: yes

Circularity Check

1 steps flagged

30% LRD duty cycle and 0.3 dex scatter are parameters tuned to reproduce observed LFs

specific steps
  1. fitted input called prediction [Abstract]
    "SMBHs with M_BH≤10 M_seed ∼10^7 M_⊙ stay in the LRD phase with a duty cycle of 30%. This simple framework reproduces the observed LFs and clustering of LRDs. ... although we find that a ∼0.3 dex log-normal scatter in bolometric luminosity is required to reproduce the bright end."

    The duty cycle fraction and scatter amplitude are stated as chosen/adjusted specifically so the model matches the observed luminosity functions; the subsequent claim that the framework 'reproduces' those LFs therefore reduces to the input choice by construction.

full rationale

The central claim that the empirical model reproduces observed LFs and clustering rests on an explicitly fitted duty cycle (set to 30% for M_BH ≤ 10^7 M_⊙) and luminosity scatter (0.3 dex) chosen so the simulated population matches the JWST LRD LF and bright-end quasar LF. This makes the LF reproduction a fitted outcome rather than an a-priori prediction from the underlying simulation. Clustering reproduction is not shown to require additional tuning and may retain some independence, but the load-bearing LF match is by construction. No other circular steps (self-citation chains or ansatz smuggling) are evident from the provided text.

Axiom & Free-Parameter Ledger

2 free parameters · 1 axioms · 0 invented entities

The central claim depends on two fitted parameters in the empirical mapping and on the assumption that the simulation's reionization physics is accurate enough to identify AGN as the dominant HeII driver.

free parameters (2)
  • LRD duty cycle = 30%
    Fraction of time SMBHs below 10^7 solar masses spend in the LRD phase; set to 30% to match observed LRD number densities.
  • bolometric luminosity scatter = 0.3 dex
    Log-normal scatter added to match the bright end of the quasar LF; set to 0.3 dex.
axioms (1)
  • domain assumption The Lumina simulation self-consistently follows ionization and thermal state of the IGM driven by the AGN population.
    Invoked when stating that the same AGN population drives HeII reionization.

pith-pipeline@v0.9.1-grok · 6018 in / 1507 out tokens · 36355 ms · 2026-06-30T15:10:09.293043+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Forward citations

Cited by 2 Pith papers

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score.

  1. Little Red Dots as Intermediate Mass, Super-Eddington Engines: Insights from Type IIn Supernovae and The 1837-1856 Great Eruption of $\eta$ Carinae

    astro-ph.GA 2026-06 unverdicted novelty 6.0

    LRDs are reinterpreted as intermediate-mass super-Eddington systems with wind-driven pseudo-photospheres that explain their spectra and imply engine masses below 10^5 solar masses rather than overmassive black holes.

  2. The Lumina Project: Intergalactic Clumping and Recombination Sinks

    astro-ph.GA 2026-06 unverdicted novelty 4.0

    Simulations show recombination-weighted clumping is systematically lower than density-based measures, density-only prescriptions overpredict rates by 1.29-1.84 depending on redshift, and a new phase-space clumping fac...

Reference graph

Works this paper leans on

1 extracted references · 1 canonical work pages · cited by 2 Pith papers · 1 internal anchor

  1. [1]

    Origins of the UV continuum and Balmer emission lines in Little Red Dots: observational validation of dense gas envelope models enshrouding the AGN

    AbramowiczM.A.,CzernyB.,LasotaJ.P.,SzuszkiewiczE.,1988,ApJ,332, 646 Agazie G., et al., 2023, ApJ, 951, L8 Aird J., et al., 2010, MNRAS, 401, 2531 Aird J., Coil A. L., Georgakakis A., Nandra K., Barro G., Pérez-González P. G., 2015, MNRAS, 451, 1892 Aird J., Coil A. L., Georgakakis A., 2018, MNRAS, 474, 1225 Akins H. B., et al., 2025a, ApJ, 980, L29 Akins ...

    work page internal anchor Pith review Pith/arXiv arXiv doi:10.5281/zenodo.1451799 1988