© The Authors 2025

1. Introduction

Over the past year, the James Webb Space Telescope (JWST) has opened a new window on the high-redshift Universe, discovering galaxies and Active Galactic Nuclei (AGNs) at unprecedented redshifts (see Adamo et al. 2024, for a review). The vast majority of these AGNs have been selected or confirmed through broad emission lines, which generally are a signpost for the deep potential well of massive black holes (MBHs, e.g., Wasleske & Baldassare 2024).

A long-standing technique to pre-select high-z AGN candidates is color selection (e.g., Fan et al. 2001). The different shape of galaxy and AGN spectral energy distributions (SEDs) modulates the colors of sources, isolating in different regions those that are AGN- or galaxy-dominated, at least when the ratio between the luminosity of the two components is sufficiently different (Volonteri et al. 2017; Valiante et al. 2018). Additionally, at optical wavelengths AGN can also be selected by variability (e.g., Hayes et al. 2024) and by having a compact appearance (e.g., Jarvis & MacAlpine 1998), as AGN appear as point sources, and when sufficiently more luminous than the host galaxy outshine it.

Search for broad-line emitters with JWST were performed on targets selected in the JADES survey (Maiolino et al. 2024b), using data from spectroscopy programs (Harikane et al. 2023; Kocevski et al. 2023) and in the ASPIRE spectroscopic survey (Lin et al. 2024). Matthee et al. (2024) performed a search for broad-line emitters in the FRESCO (Oesch et al. 2023) and EIGER (Kashino et al. 2023) finding an abundant population, characterized by red colors and compact sizes, winning them the name “little red dots”. Greene et al. (2024) searched in the UNCOVER survey (Bezanson et al. 2024) for compact red sources, and follow-up spectroscopy confirmed that they were broad line emitters. Many if not most of the little red dots show a blue excess in the short-wavelength NIRCAM bands (Greene et al. 2024; Kocevski et al. 2024), which has been ascribed to either the host galaxy emission or scattered light from the AGN. Following these discoveries, other photometric searches were performed, selecting for red sources (e.g., Kokorev et al. 2024; Kocevski et al. 2024) and identifying a very large population of sources. If all these red sources are AGN, they would severely challenge our understanding of galaxy and massive black hole (MBH) evolution (Akins et al. 2024).

In this paper, we want to investigate the properties of MBHs in galaxies to understand the relationship between little red dots and the general population of AGN, as well as what is needed for an AGN to be a little red dot. Specifically, we model and analyze simulated photometry of galaxies and AGN based on the OBELISK simulation, select red sources, which we refer to as “color-selected little red dots” (cLRDs), and explore their properties.

2. Methodology

2.1. The OBELISK simulation

OBELISK (Trebitsch et al. 2021) is a high-resolution resimulation down to redshift z ≃ 3.5 of the most massive halo at z = 1.97 in the HORIZON-AGN (Dubois et al. 2014) simulation. The simulation is based on a Λ CDM cosmology with WMAP-7 parameters (Komatsu et al. 2011): H₀ = 70.4 kms⁻¹ Mpc⁻¹, Ω_Λ = 0.728, Ω_m = 0.272, Ω_baryon = 0.0455, σ₈ = 0.81, and spectral index n_s = 0.967. The zoomed-in region has a volume of ∼10⁴ h⁻³ cMpc³. We summarize here the main assumptions, and refer the readers to Trebitsch et al. (2021) and Dong-Páez et al. (2023b) for further details.

OBELISK was run with RAMSES-RT (Rosdahl et al. 2013; Rosdahl & Teyssier 2015), a radiative transfer hydrodynamical code that builds on the adaptive mesh refinement (AMR) RAMSES code (Teyssier 2002). A cell is refined up to a smallest size of 35 pc if its mass exceeds eight times the mass resolution. The simulation assumes an ideal monoatomic gas with an adiabatic index of γ = 5/3 and includes gas cooling and heating down to very low temperatures (50 K) with nonequilibrium thermochemistry for hydrogen and helium, and a contribution to cooling from metals (at equilibrium with a standard ultraviolet background) released by supernovae (SNe).

In the high-resolution region, dark matter particles have masses of 1.2 × 10⁶ M_⊙, star particles have masses of ∼10⁴ M_⊙, and assume a Kroupa (2001) initial mass function between 0.1 and 100 M_⊙. Stars form in gas cells with a density higher than 5 H cm⁻³ and Mach number of ℳ ≥ 2, with a star formation efficiency that depends on the local gas density, sound speed, and turbulent velocity. Mechanical SN feedback is released 5 Myr after the birth of a stellar particle, with a mass fraction of 0.2, and releasing 10⁵¹ erg per SN. OBELISK also includes modeling of dust as a passive variable.

Massive black holes form when in cells where both gas and stars exceed a density threshold of 100 H cm⁻³ and the gas is Jeans-unstable, with an exclusion radius of 50 comoving kpc to avoid the formation of multiple MBHs. The initial MBH mass is 3 × 10⁴ M_⊙, and MBHs grow by accretion, capped at the Eddington rate, using a Bondi-Hoyle-Lyttleton (BHL) formalism in which the accretion rate is calculated via the local average gas density, gas sound speed, and MBH relative velocity with respect to the gas.

OBELISK includes a dual-mode AGN feedback. Below f_Edd <  0.01 (where f_Edd is the ratio of the MBH mass accretion rate over its Eddington mass accretion rate), an AGN releases a fractio,n ε_MCAF, of the accretion energy as kinetic energy in jets. Jets assume magnetically chocked accretion flows, and ε_MCAF is a polynomial fit to the simulations of McKinney et al. (2012) as a function of MBH spin. At higher f_Edd ≥ 0.01, the 15% of the accretion luminosity is released isotropically as thermal energy. MBH spin magnitudes and directions are self-consistently evolved on the fly via gas accretion and MBH-MBH mergers. Two MBHs are merged when their separation becomes smaller than 4Δx. The simulation models dynamical friction explicitly, including both gas and collisionless particles (stars and DM) (following the implementation presented in Dubois et al. 2013; Pfister et al. 2019). The dark matter distribution is smoothed by a cloud-in-cell interpolation to avoid artificial scattering.

OBELISK includes a model for dust evolution. The model assumes that dust grains are released in the ISM via SN explosions, grow in mass via accretion of gas-phase metals, and are destroyed by SN explosions and via thermal sputtering. Dust grains are further assumed to be perfectly coupled to the gas, so that they can be described by a scalar passively advected like the metallicity. We assume that all grains have an average size of a_g = 0.1 μm and density μ_g = 2.4_g cm⁻³. Further details can be found in Trebitsch et al. (2021) and Dubois et al. (2024).

The galaxy and AGN populations in OBELISK are presented in Trebitsch et al. (2021). Galaxies are consistent with the expected relation between stellar and dark matter halo mass, and the star formation rate density is in good agreement with that found in overdense regions. Galaxy masses are also not challenged by JWST results (Keller et al. 2023). The simulation broadly produces gas, metals, and dust content in high-z galaxies in agreement with observed scaling relations with stellar mass (Trebitsch et al. in prep.). The properties of the galaxies analyzed in this work are shown in Fig. 1. AGN in OBELISK are abundant but MBHs mainly grow only in galaxies with M_star >  10⁹ M_⊙ because of the effect of SN feedback, and in part also because of AGN feedback limiting their own growth. At z = 6 the most massive MBH has M_BH = 1.3 × 10⁶ M_⊙. As a consequence, OBELISK’s AGN do not reach very high luminosity, the maximum luminosity being 3.5 × 10⁴³ erg s⁻¹.

Fig. 1. Properties of simulated galaxies with Mstar >  108 M⊙ at z = 6, including galaxies that do not host MBHs. From top to bottom: star formation rate averaged over 10 Myrs in black, and over 100 Myrs in gray; mass-weighted stellar age, and mass-weighted stellar metallicity, where solar metallicity is assumed to be 0.02, a dust mass within twice the effective radius. Reff is defined as the half-mass projected radius in the plane orthogonal to the stellar angular momentum, computed as the (2D) radius such that 50% of the mass is within that radius. The violet line in the top panel shows the fit from Speagle et al. (2014).

2.2. Galaxy and AGN absorption and photometry

A sketch of our approach is shown in Fig. 2. Besides intrinsic galaxy properties, such as, stellar mass (M_star), SFR, dust-to-gas ratio, photometry of the galaxies in OBELISK has also been calculated in NIRCAM bands (F115W, F150W, F200W, F277W, F356W, F444W), both considering the intrinsic emission determined by the star formation history and by attenuating it, as described in Trebitsch et al. (in prep.). Stellar population SEDs are based on Bpass v2.2.1 binary stellar population model (Eldridge et al. 2017; Stanway & Eldridge 2018). They depend on the age and metallicity of the stellar population. The intrinsic SED is sampled using 10⁶ photons that probe the distribution of the luminosities of, and are assigned to, star particles, and each photon is assigned a restframe wavelength λ in the range 0.075−1 μm based on the SED of its emitting star particle, as well as an isotropically selected initial direction. Photons are then propagated using the RASCAS code (Michel-Dansac et al. 2020), based on the optical depth along the line of sight, computed using the dust mass density in each cell. In the case where a photon is absorbed, then the propagation stops, otherwise the photon is scattered in a new direction. Photons are collected in 12 independent directions to construct the observed SEDs.

Fig. 2. Sketch of our approach to modeling gas and dust absorption, showing how differential attenuation affects different parts of a given galaxy and lines of sight. We highlight two lines of sight with rays shown in purple and red arrows for each emitter. Darker gray colors show gas and dust clouds with larger optical depth. Stars’ colors indicate age: for instance young (blue) stars are generally still embedded in their birth clouds, while this is not necessarily the case for older (yellow and orange) stars. Since at least some star-forming regions are dusty, their contribution to the SED will be more affected by the extinction curve than other parts of the galaxy. The position of the AGN is marked in black, with thicker arrows. The example distribution of emitters and absorbers shown here has been chosen to highlight how the same source can be AGN- or stellar-dominated, depending on the line of sight. In this example, along the cumulative purple line of sight the AGN dominates over starlight: the AGN is seen through a clear line of sight, while much of the emission from the stellar population is absorbed. Vice-versa, along the cumulative red line of sight the AGN emission is sub-dominant to starlight.

A full discussion of galaxy dust properties and comparison with observations is presented in Trebitsch et al. (in prep.). We here note only that the dust masses of z = 4 − 9 in OBELISK reproduce well the observed scaling and normalization with stellar mass. In terms of dust attenuation, the typical A_V increases from about 0.3 to 1.5 − 2.0 for galaxies with masses from 10⁸ M_⊙ to 10¹⁰ M_⊙ (depending on the choice of the extinction curve), with a dispersion of about 0.7 mag around this value. Finally, the observed UV β slopes (not accounting for nebular emission) range from β ≃ −3 to −0.5 assuming a SMC-like extinction curve, and to −1.5 for a MW-like extinction curve. We note that even for a given extinction curve, OBELISK reproduces a variety of attenuation curves due to the complex geometrical effects between the stellar and dust distributions.

The AGN continuum is described by the following equation (Volonteri et al. 2017):

$$\begin{array}{c}\hfill f_{\nu }=\mathcal{N}({\nu }^{{\alpha }_{\mathrm{UV}}}{e}^{-\frac{h\nu }{k{T}_{\mathrm{peak}}}}{e}^{-\frac{k{T}_{\mathrm{IR}}}{h\nu }}+a{\nu }^{{\alpha }_{\mathrm{X}}}),\end{array}$$(1)

with ν the photon frequency, h the Planck constant, α_UV = −0.5, α_X = −1, and kT_IR = 0.01 Ryd. The continuum is based in optical/near-IR on the Shakura-Sunyaev solution (Shakura & Sunyaev 1973), following Thomas et al. (2016) for the estimation of T_peak. T_peak represents the location of the peak of the spectrum, and it occurs at temperature of the order of 3 × 10⁶ − 10⁶ − 7 × 10⁴ K for black holes with f_Edd = 0.1 and mass 10⁴ − 10⁶ − 10⁹, respectively. As per the standard Shakura-Sunyaev solution, T_peak increases with decreasing black hole mass and increasing f_Edd. The formalism in Thomas et al. (2016) allows one to obtain a good match with standard bolometric corrections in the luminosity limit where the corrections have been derived. We refer the reader to Fig. 1 and Appendix A in Volonteri et al. (2017) for additional information. The model is calibrated in X-rays using results from the physical models developed by Done et al. (2012): the normalization a is obtained through α_OX, the exponent of a power-law connecting the continuum between 2 keV and 2500 Å, fitting the dependence on MBH mass and Eddington ratio using the results in Dong et al. (2012). The last term in Equation (1) describes the corona contribution and it is set to zero below 0.1T_peak, where the contribution becomes negligible: this avoids that the power-law pops up at low frequencies below the peak at T_peak. The global normalization 𝕉 is obtained by requiring that the bolometric luminosity matches the accretion luminosity L_bol = 1.26 × 10³⁸ f_Edd(M_BH/M_⊙) erg s⁻¹. X-ray luminosity is calculated in the [2−10] keV range (observer frame). We have not attempted to include scattered light contributing to the UV emission.

We consider two extinction laws: one similar to the Milky Way (“MW-like”) with 54% of the grain number density in silicate grains and 46% in carbonaceous grains, and one similar to the Small Magellanic Cloud (“SMC-like”) with 100% silicate grains. They both use the Mathis et al. (1977) grain size distribution ranging from 0.001 μm to 0.25 μm (slope −3.5), and dust grain optical constants from Weingartner & Draine (2001) and Laor & Draine (1993). For a given amount of dust, the MW-like extinction curve, compared to the SMC-like, leads to more extinction overall but to much less reddening (i.e., to less extinction in the blue part, λ ≲ 0.2 μm restframe, relative to the red part of the spectrum, λ ≳ 0.25 μm restframe). To model absorption for AGN we calculate the contribution of the interstellar medium (ISM) to the column density at the location of the MBHs by casting rays isotropically around each MBH in the outputs of the simulation. We used a version of the RASCAS code (Michel-Dansac et al. 2020) modified to integrate column densities very efficiently. Values of 12 individual rays per galaxy, for both the stellar population and the AGN, are employed throughout the rest of the paper unless stated otherwise. To estimate the AGN attenuation at JWST wavelengths we multiply the hydrogen column density by the mean dust-to-gas mass ratio within 2R₅₀, where R₅₀ is the 3D half-mass radius, to obtain the dust column density and then convolve it with the extinction laws to obtain the optical depth. We have checked that there is no significant gradient in the dust-to-gas mass ratio between R₅₀ and 2R₅₀. The typical dust-to-gas mass ratios are between 10⁻² and 10⁻³.

To estimate the AGN gas column density, we first integrated the gas column density from the MBH position to the virial radius of its host halo for each of the 12 rays (‘max’ column density). This results in very high gas column densities, making the AGN contribution in JWST bands completely negligible and producing a very high fraction of obscured AGN in X-rays. Indeed, the gas distribution near the MBHs is stochastic – depends on the position of few gas cells – and the inclusion of boosted gas dynamical friction makes the MBH cling to dense gas. Furthermore, simulations find that the column density around an AGN decreases after feedback has been able to clear out paths in its surroundings (Trebitsch et al. 2019; Ni et al. 2020; Lupi et al. 2022). In OBELISK MBHs start to grow when galaxy masses reach M_star ∼ 10⁹ M_⊙, this is mainly because SN feedback makes gas very turbulent in low-mass galaxies. This retarded growth gave a better match to pre-JWST AGN luminosity functions (Trebitsch et al. 2021, see also Fig. 4), but as MBHs struggle to grow they have not been able to open these chimneys in most of OBELISK’s galaxies. To account for these considerations, we base our fiducial “outmax80” model on the proposal of Fabian et al. (2008) that the radiative pressure on dust evacuates dusty gas in the MBH vicinity at high Eddington ratios, or that it prevents its accumulation (Fabian et al. 2009). Fabian et al. (2008) show that for a given Eddington ratio there is a minimum column density for gas to counteract radiation pressure and not be blown away. In practice, this provides a limiting column density value as a function of the Eddington ratio below which dusty gas will be removed or prevented to accumulate. We include this effect in the model as follows. We fit the curve of this threshold column density as a function Eddington ratio (Fabian et al. 2008) and compare it to the column density seen by the MBH. If the “max” column density is below this threshold value, we consider that the inner region has been evacuated and integrate the column density outside a sphere of radius¹ 80 pc. Otherwise we retain the full “max” column density. In the following we consider this as the fiducial case for the column density, since it captures the effect that feedback from actively growing MBHs would have on the galaxy, effect that is not present in OBELISK’s galaxies because of the limited MBH growth.

To estimate the X-ray obscuration properties, we use the hydrogen column density, and we furthermore consider that dust-free gas can linger between the edge of the accretion disk, calculated as the self-gravity radius, and the dust sublimation radius. As noted by Maiolino et al. (2024a) the mismatch between high hydrogen column densities and low dust column densities is a typical feature even in low-redshift AGNs, and the addition of this gas eases the tension between the X-ray and optical AGN populations. However, one has to justify that just outside the dust sublimation radius radiation does not find a wall of dusty gas, which would happen if the density distribution were continuous. The “outmax80” model we adopt, where radiative pressure on dust evacuates dusty gas, would explain how a very high gas density within the dust sublimation radius can drop to low dusty gas densities immediately afterward. To include the unresolved contribution below the dust sublimation radius we follow the approach of Dong-Páez et al. (2023a). Assuming that gas falls radially under the effect of the MBH gravity and a spherically symmetric configuration the gas density in the torus should follow a power-law density profile ρ ∝ r^−3/2. The density profile is normalized so that the total mass between the Bondi radius and the edge of the accretion disk, which we assume corresponds to the self-gravity radius, corresponds to the accretion rate multiplied by the free fall time. For the self-gravity radius we adopt the expression from (Laor & Netzer 1989):

$$\begin{array}{c}\hfill r_{\mathrm{sg}}\simeq 3586r_{\mathrm{g}}{\alpha }^{2/9}{\left(\frac{M_{\mathrm{BH}}}{{10}^8{\mathrm{M}}_{\odot }}\right)}^{-2/9}{f}_{\mathrm{Edd}}^{4/9},\end{array}$$(2)

where r_g = GM_BH/c², G is the gravitational constant, c is the speed of light in vacuum, and we set α = 0.1. We then integrate the density profile from r_sg to the dust sublimation radius (Suganuma et al. 2006):

$$\begin{array}{c}\hfill r_{\mathrm{subl}}\approx 0.47{\left(\frac{L_{\mathrm{UV}}}{{10}^{46}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{1/2}\mathrm{pc}.\end{array}$$(3)

The final expression is:

$$\begin{array}{c}\hfill N_{\mathrm{H},\mathrm{inner}}=\frac{3\sqrt{2}\pi }{4{\epsilon }_{r}c{\sigma }_{\mathrm{T}}}{f}_{\mathrm{Edd}}{(G{M}_{\mathrm{BH}})}^{1/2}{r}_{\mathrm{subl}}^{-1/2}[{\left(\frac{r_{\mathrm{sg}}}{r_{\mathrm{subl}}}\right)}^{-1/2}-1],\end{array}$$(4)

where ε_r is the radiative efficiency and σ_T is the Thomson cross section. For ease of implementation we fit N_H, inner from the sample of MBHs in galaxies with mass M_star >  10⁸ M_⊙ and f_Edd >  10⁻⁴ at z = 6 to obtain:

$$\begin{array}{c}\hfill log\left(\frac{N_{\mathrm{H},\mathrm{inner}}}{{\mathrm{cm}}^{-2}}\right)=0.72log(f_{\mathrm{Edd}})-0.29log\left(\frac{M_{\mathrm{BH}}}{{10}^8{\mathrm{M}}_{\odot }}\right)+24.90,\end{array}$$(5)

with RMS residual = 0.16 when we consider the actual value of ε_r from OBELISK, which includes spin evolution, and in which spin increases with MBH mass (so radiative efficiency decreases), as discussed in Dong-Páez et al. (2023b). If we consider a fixed ε_r = 0.1, we find

$$\begin{array}{c}\hfill log\left(\frac{N_{\mathrm{H},\mathrm{inner}}}{{\mathrm{cm}}^{-2}}\right)=0.77log(f_{\mathrm{Edd}})+0.10log\left(\frac{M_{\mathrm{BH}}}{{10}^8{\mathrm{M}}_{\odot }}\right)+23.13,\end{array}$$(6)

with RMS residual = 0.001, but in this case all AGN with L_{2 − 10keV} >  3 × 10⁴³ erg s⁻¹ would have N_H, inner above 10²³ cm². We note, however, that our model of this dust-free gas is very simplistic, assuming spherical symmetry and a density power-law slope of −3/2, therefore we consider this more of an indication of the possible magnitude of the effect than an actual quantitative estimate. For the remainder of the paper we consider Eq. (5) when we include unresolved obscuration.

2.3. Massive black hole population

The “true” MBH properties in the simulation do not contain any overmassive MBHs that can be compared to JWST AGN candidates, and many MBHs accrete at low Eddington ratios, therefore we create a population, “distr” where the MBH masses are extracted from the M_BH-M_star relation from Greene et al. (2020, G20):

$$\begin{array}{c}\hfill log\left(\frac{M_{\mathrm{BH}}}{{\mathrm{M}}_{\odot }}\right)=6.87+1.39log\left(\frac{M_{\mathrm{star}}}{{10}^{10}{\mathrm{M}}_{\odot }}\right).\end{array}$$(7)

We also consider an intrinsically overmassive population (’distr overm’) based on the fit by Pacucci et al. (2023, P23):

$$\begin{array}{c}\hfill log\left(\frac{M_{\mathrm{BH}}}{{\mathrm{M}}_{\odot }}\right)=8.17+1.06log\left(\frac{M_{\mathrm{star}}}{{10}^{10}{\mathrm{M}}_{\odot }}\right).\end{array}$$(8)

By mocking the MBH population to obtain sufficiently luminous AGN we can address what type of MBH properties are needed to have cLRDs, given a normal population of galaxies. We cannot, however, address how cLRDs are born: that would require self-consistent modeling linking MBH and galaxy evolution, in particular the star formation history and how AGN emission and outflows modulate the dust properties (e.g., Elvis et al. 2002; Gaskell et al. 2004; Antonucci 2023). In practice, all properties related to galaxies and the interstellar medium are self-consistent, while for the AGN we retain only its position in the galaxy, which determines the gas and dust column density seen by the mocked AGN.

Based on Volonteri et al. (2017) and Greene et al. (2024), Eddington ratios are extracted from a normal distribution centered in log(0.25) with standard deviation 0.5 dex for “distr” and “overm”, and centered in log(0.5) for a “sEdd” case, where if f_Edd >  0.3 we shift redward the energy peak of the spectrum to mimic photon trapping, using r_cor = 20 for sub-Eddington and r_cor = 200 for super-Eddington (see Kubota & Done 2019, for a fully physical model). Here r_cor is in units of the MBH gravitational radius, and sets the peak temperature of the SED. Its value in the sub-Eddington case has been calibrated to match standard bolometric corrections in the relevant limits (Volonteri et al. 2017). In either case we do not cap f_Edd and we do not decrease the overall luminosity for sources with f_Edd >  0.3, but we have checked that limiting the luminosity to L_Edd does not change quantitatively our results. The effect that quantitatively changes the results in the sEdd case is related to the relative increase in luminosity of the AGN with respect to the galaxy at the redder JWST bands (356, 444) with respect to the bluer ones (115, 200) because the energy peak in the sEdd case shifts toward lower frequencies.

Model assumptions are summarized in Table 1. Since this analysis is mainly focused on AGN with bolometric luminosity > 10⁴⁴ erg s⁻¹ we consider galaxies in the OBELISK simulation at redshift z = 6 and M_star >  10⁸M_⊙. Although we post-process M_BH and f_Edd we analyze galaxies hosting MBHs in the simulation, unless specified otherwise. The reason is that we need to know the MBH position in order to calculate attenuation. For the analysis of the absorbed fraction we extend the analysis to z = 4 and z = 5. For the bolometric luminosity function, we consider all galaxies with M_star >  10⁸ M_⊙ at z = 5 and z = 7. In future work we will explore the redshift evolution of AGN and cLRDs.

Finally, at z = 6 the simulation includes 340, 75 and 8 galaxies with M_star = 10⁸ − 10⁹ M_⊙, 10⁹ − 10¹⁰ M_⊙ and > 10¹⁰ M_⊙, respectively. This means that the high-mass end of the sample is sparsely populated, and the results will be highly dependent on the randomly picked M_BH and f_Edd for the mock sample: drawing a single M_BH and f_Edd per galaxy would not sample the full parameter space. We therefore draw a different M_BH and f_Edd from the respective distributions for each of the 12 lines of sight in order to obtain a sufficiently large sample, except in Sect. 6, where we explore the effect of line-of-sight differences. At z = 4 and z = 5, used only to estimate the X-ray absorbed fraction, we draw 2 and 4 cases per galaxy, using the mean column density rather than 12 lines of sight, in order to have a similar sample size at each redshift.

To validate the model and show that it produces a reasonable AGN population, the AGN bolometric luminosity functions at z = 5 and z = 7 (for M_star >  10⁶ M_⊙ to avoid incompleteness) are shown in Fig. 4. To obtain the cosmic average we have rescaled the OBELISK galaxy mass function, which simulated an overdensity, to that of the NEWHORIZON simulation, which simulated an average region of the Universe. The rescaling of the mass functions provides a correction for the number density of galaxies as a function of their mass and redshift. To derive the rescaled AGN LF, we calculate the individual AGN LFs in different M_star bins, renormalize the number density in each mass bin based on the correction calculated from the ratio of the NEWHORIZON and OBELISK galaxy mass function in that mass bin, and then sum the individual AGN LFs obtained in the different M_star bins to obtain the complete LF. We have checked that we obtain similar results rescaling by the halo mass function, and binning in halo mass, but the number density correction has a steeper dependence on halo mass and therefore we opted for rescaling based on galaxy mass to have a smoother evolution.

The simulated AGN LF is compared in Fig. 4 to the bolometric luminosity function of AGN from Matthee et al. (2024, M+24), Greene et al. (2024, G+24) and Shen et al. (2020, S+20). The observational estimate of the luminosity function (LF) around becomes more uncertain as the luminosity decreases. The spectroscopic analysis in FRESCO and EIGER in M+24 did not detect any broad line AGN with bolometric luminosity larger than ∼10⁴⁶ erg s⁻¹, therefore this gives an upper limit at high luminosity (yellow downward arrows in Fig. 4). On the other hand, these bolometric LFs have been derived from broad line AGN only, and Scholtz et al. (2023) finds a comparable-size population of narrow-line AGN. The extrapolations of some of the LFs would likely overpredict the bright end of the LF as known pre-JWST. We note, however, that for ease of interpretation in Fig. 4 and throughout the paper we have assumed an active fraction of unity: each galaxy hosts an MBH that is accreting at a non-negligible Eddington ratio. This is therefore an upper limit. Assuming that only 10% of MBHs independently of MBH or galaxy mass are active, for instance, would decrease the LF by a factor of 10. We discuss how an active fraction of 100% translates into different detectable AGN fractions in Sect. 5.

3. Spectral energy distributions and absorption

At these redshifts, for our unattenuated models the galaxy dominates the emission at short JWST wavelengths, and the AGN at long JWST wavelengths (see Volonteri et al. 2017, 2023, for models without dust), as long as the AGN is sufficiently bright with respect to the host galaxy. This can be intuitively understood as follows (see Fig. 2 in Volonteri et al. 2017). In absolute terms, galaxies are star-forming and dominated by young stellar populations (Fig. 1), with a blue intrinsic spectrum that falls off more steeply than the AGN intrinsic SED at longer wavelengths. Furthermore, AGN SEDs peak at lower frequencies as the MBH mass decreases, and the AGN emission dominates only when the MBH represents a sizeable fraction of the galaxy mass. The dominance of galaxy emission over the AGN at short JWST wavelengths is only accentuated when including attenuation.This is shown in Fig. 3 for some idealized cases. There is, however, a variety of behaviors, and we will show more realistic examples in Sect. 6. In the left panels of Fig. 3 we vary smoothly M_BH and f_Edd at fixed gas and dust column density seen by the AGN, and show random line of sight through the galaxy for the stellar population. In the right panels we fix M_BH and f_Edd, vary smoothly the AGN column density, and show three different lines of sight for the stellar population. We can see that in this model the wavelength where the SEDs transition from the galaxy to AGN dominated changes as a function of dust model, column density and MBH properties. Those that are selected as cLRDs have the transition in the “right place” to obtain a v-shaped SED.

Fig. 3. Examples showing how different assumptions affect the emergent SEDs at z = 6. For the same galaxy we show different MBH, fEdd and NH, as well as different lines-of-sight for the emission from the stellar population. The NIRCAM response for the filters analyzed in this paper is shown in gray.

Fig. 4. Bolometric luminosity function for the actual MBHs in the simulation (“true”) and for the models described in Table 1, where “resc” means we have corrected the number densities for the overdensity bias: OBELISK is a protocluster. To obtain the cosmic average we have rescaled the OBELISK galaxy mass function to that of the NEWHORIZON simulation, which simulates an average region of the Universe. This figure assumes an active fraction of unity, and therefore represents an upper limit to the LF for each of the models. The results are compared to the luminosity functions derived in Matthee et al. (2024, M+24) and Greene et al. (2024, G+24) from JWST data and to fits to the evolution of the bolometric luminosity function derived in Shen et al. (2020, S+20, pre-JWST).

Fig. 5 shows the dependence of the mean hydrogen column density from the ISM only on galaxy mass and size. In general low-mass galaxies are diffuse because of efficient SN feedback puffing up the gas component (hence stars follow). While at larger masses (M_star ≳ 10⁹ M_⊙), SNe become inefficient, gas can pile up, and leads to wet compaction (Dekel & Burkert 2014). We first discuss galaxy sizes: 75% of galaxies with M_star >  10⁹ M_⊙ have half mass radius smaller than 300 pc (100% for M_star >  3 × 10⁹ M_⊙). The most compact galaxy in our sample has R_50, mass = 55 pc for a mass of M_star = 1.84 × 10⁹ M_⊙. The “max” column density is clearly highly correlated to galaxy mass and compactness, while, obviously, removing a sphere with a radius of 80 pc makes the most compact galaxies less optically thick. This is particularly noticeable for galaxies with R_50, mass <  100 pc. In Fig. 5 one can also appreciate the importance of the inhomogeneites in gas distribution: for a subset of galaxies we show the 12 lines of sight we use, and we see that for a given galaxy there can be a spread of 1-2 orders of magnitude. The ISM column density seen by the AGN increases with stellar mass as galaxies become more compact. A contribution of 10²² − 10²³ cm⁻² from the ISM to the total column density is in good agreement with the results obtained by Silverman et al. (2023) on 3 AGNs at z ∼ 0.75 − 2.

Fig. 5. Relation between hydrogen column density in the ISM and galaxy mass and half-mass radius at z = 6. Top panel: column density obtained integrating the column density from the position of the MBH to the host galaxy virial radius (“max” model). Bottom panel: column density obtained excluding the inner 80 pc around the MBH (“outmax80” model). The orange points show the mean column density for each galaxy, with size scaling with Mstar. For a random subset of 24 galaxies we show the column density for 12 different lines of sight, with the grayscale proportional to Mstar, as shown in the color bar.

In the bottom panels of Fig. 6 we show “outmax80” AGN gas column densities for models that produce more than 10% of cLRDs (see Sect. 4). The bifurcation at M_star > 10¹⁰ M_⊙ is related to the ability of radiation pressure to lift gas: the high branch is where the column density is too large for radiation pressure on dust to lift dusty gas, the low branch is the opposite case. This is noticeable only at the high mass end because in most low mass galaxies the intrinsic column density is below the limit to counteract radiation pressure. In general, most of our simulated AGN have gas column density above 10²²cm s⁻², therefore unobscured, UV-bright AGN are a minority, associated with favorable lines of sight. Simulations typically find low average column densities only after the energy injection from the AGN has cleared out its surroundings (Trebitsch et al. 2019; Ni et al. 2020; Lupi et al. 2022). However, since in Obelisk the MBHs have not grown much, they have injected little energy.

Fig. 6. Hydrogen column density as a function of galaxy stellar mass. Top three rows: distribution of “effective” hydrogen column density for starlight obtained from the ratio of intrinsic and attenuated fluxes in each JWST filter we have considered (Eq. (9)). Bottom row: distribution of hydrogen column density seen by the AGN as a function of galaxy stellar mass at z = 6. This is the column density integrated above the dust sublimation radius, for example the column density that enters into calculating dust absorption. The red points highlight the cases that host a cLRD meeting the UNCOVER selection.

The galaxy SED is modulated by the star formation history of the galaxy and by the gas and dust distribution. Attenuation for starlight is not the same at all wavelengths, but it depends on where stars are located: for instance, young stars dominating the emission at short wavelengths are typically embedded in dense gas clouds, and therefore they are surrounded by higher column densities. This can be appreciated in the top panels of Fig. 6, where we have rescaled the optical depth by the cross section and dust-to-gas ratio to obtain an effective column density:

$$\begin{array}{c}\hfill N_{\mathrm{H}}=\frac{X_{\mathrm{H}}{\tau }_{\lambda }}{{\kappa }_{\lambda }{m}_{\mathrm{p}}\mathrm{DTG}},\end{array}$$(9)

where τ_λ is obtained from the intrinsic and attenuated fluxes at wavelength λ, κ_λ is the absorption coefficient, m_p is the proton mass, X_H = 0.76 is the universal hydrogen mass fraction and DTG is the dust-to-gas mass ratio. This effective column density can be compared to the AGN’s, shown in the bottom panels. The AGN is generally much more absorbed than starlight, since by definition it is located at the core of the galaxy, and this is the case for cLRDs as well, as shown by the red points in all panels. This also explains why the dust model, MW versus SMC, has a stronger impact on the AGN compared to the galaxy. The ratio of attenuated luminosities for two dust models at fixed intrinsic luminosity is L_att, dust1/L_att, dust2 = exp(N_HDTG(κ_dust2−κ_dust1)m_p/X_H). Since N_H is larger for AGN than for starlight, the difference between two dust models is stronger for AGN than for the stellar population.

The X-ray absorbed fraction, defined as the fraction with N_H >  10²³ cm⁻², is shown in Fig. 7. Except for the “max” model, which overpredicts the absorbed fraction, all other models underpredict the absorbed fraction from X-ray observations, unless we include the hydrogen column density within the sublimation radius. We have not attempted to modify the X-ray emission to account for the lack of X-ray emission from some/many high-z AGN and red quasars (Ma et al. 2024a; Yue et al. 2024; Ananna et al. 2024); 2024arXiv240500504M. Given that we find hydrogen column densities below 10²⁵ cm⁻² the X-ray luminosity of the AGN at 2-10 keV observer’s frame, corresponding to ∼10 − 40 keV rest frame, would be dimmed by at most by 10% for L_{2 − 10keV} >  10⁴³ erg s⁻¹; this is small compared to the discrepancy seen in observations that is more than 2 orders of magnitudes for most JWST AGN. In the L_{2 − 10keV} >  10⁴³ erg s⁻¹ luminosity range the typical total hydrogen column density is between 10^22.5 and 10^23.5 cm⁻². The luminosities of the simulated AGN reach 10⁴⁵ erg s⁻¹, therefore the AGN would be detectable in deep X-ray observations.

Fig. 7. Absorbed fraction, defined as the fraction of AGN with NH >  1023 cm−2 for various models compared to Vito et al. (2018) and Pouliasis et al. (2024). We include the 68% confidence region as shaded areas. The top panel shows the results including only the ISM gas, the bottom panel adds gas within the dust sublimation radius. For the cases including gas within the dust sublimation radius, we show the case with εr = 0.1, including spin evolution the absorbed fraction decreases to 60–70% (distr and sEdd), and 40% (overm).

To explain the lack of X-rays from high-z AGN one option is that X-ray and optical observations select intrinsically different populations. In this case, the X-ray bright population would be missed in optical searches because the AGN are too much reddened. However, the larger number density of broad-line AGN compared to that of X-ray selected AGN implies that the large majority of AGN is intrinsically X-ray weak (Ma et al. 2024a), or that absorption is dominated by dust-free gas, as has been suggested by Maiolino et al. (2024a). In our model, however, we do not reach column densities within the sublimation radius that would dim the X-rays at the rest frame energies: for ∼10 − 40 keV rest frame a column density > 10²⁵ cm⁻² would be needed. Another possibility is that some AGN (but not all, see e.g., Suh et al. 2024; Maithil et al. 2024) in the super-Eddington regime are intrinsically X-ray weak, because of suppression of the production of UV photons that can be upscattered due to radiation trapping, or disk winds shielding or disrupting the corona or failed thermalization (Proga 2005; Begelman & Volonteri 2017; Dai et al. 2018; Jiang et al. 2019; Pacucci & Narayan 2024). This ties in particularly well with the “sEdd” case, and is in good agreement with the suggestion that many if not most of the high-z AGN detected by JWST are in the super-Eddington regime (Schneider et al. 2023; Lupi et al. 2024; Juodžbalis et al. 2024).

4. Photometric selection of little red dots

Our model does not include emission lines, therefore we start by comparing it with photometrically selected AGN. Specifically, we apply the color selection adopted in the UNCOVER² survey (Labbe et al. 2025; Greene et al. 2024), plus a threshold in bolometric luminosity at 10⁴⁴ erg s⁻¹ to mimic the Hα flux/EW requirement. We do not impose a compactness criterion. For simplicity, we approximate the flux limit as m444W < 28, which is applied to the combined AGN+galaxy flux unless noted otherwise. The color selection requires the following criteria: −0.5 <  m115W − m200W <  1.0, m277W − m444W >  1.0, m277W − m356W >  0.7. We define as full AGN population systems with combined AGN+galaxy m444W < 28 and bolometric luminosity > 10⁴⁴ erg s⁻¹.

We first discuss the galaxy properties, without including AGN, shown in Fig. 8 (top and middle panels). The full galaxy population is shown with black points, galaxies brighter than m444W=28 in gray. None of the “galaxy-only” models get close to the cLRD classification, therefore if the cLRDs identified with this color selection are actually not AGNs, we still cannot explain their SED with our models. The same is true for models with no dust, with or without AGN (bottom panels). Certainly, our model lacks emission lines and therefore galaxies that can be scattered into the selection thanks to the presence of strong emission lines. This is however not the most likely reason. OBELISK galaxies are for the most part on the galaxy main sequence (Fig. 1, top) and do not have extreme ages, either very young or very old. The simulation therefore does not include very evolved stellar populations as in de Graaff et al. (2024), or strong starbursts (Pérez-González et al. 2024; Akins et al. 2024). If the number density of non-AGN little red dots is about 10⁻⁴ Mpc⁻³ (Akins et al. 2024), we should find 𝖪(1) of these galaxies in the OBELISK volume, therefore it’s not greatly surprising that we we do not find any. Finally, we have also assumed standard dust attenuation laws.

Fig. 8. Colors for galaxies only for SMC- (top) and MW-like (middle) dust extinction curves, and for galaxies and AGN without including dust (bottom). Galaxies at z = 6 do not enter the UNCOVER color selection highlighted by the shaded red area, nor do AGN without inclusion of dust. We consider here the “overm” case, which leads to the largest AGN contribution to colors.

Moving on to the properties of AGN selected as cLRDs, in order to skew sufficiently the colors toward red, the AGN must not contribute much (or at all) at F115W and F200W, and should instead contribute, or even dominate, at F356W and especially at F444W. This can be achieved only if the AGN is intrinsically sufficiently luminous with respect to the host galaxy: it requires MBHs to be overmassive and/or have high Eddington ratios as noted in Fig. 9. As shown in Fig. 8, using the UNCOVER selection, no AGN in any of our models, including “overm” and “sEdd” enters the selection if dust is not included. When the AGN is too luminous with respect to the host galaxy it falls out of the UNCOVER selection, though, because it becomes too red in m115W-m200W, especially in the case with SMC-like dust. The presence of emission lines, which are not included in our model, can temper these results, as they can boost the brightness in some JWST bands more than in others, for instance Hα at z ∼ 5 can boost flux in F356W and/or F444W, and help enter the selection.

Fig. 9. Colors for SMC-like dust. AGN are defined as sources having Lbol >  1044 erg s−1 and m444W < 28 (combined galaxy+AGN). cLRDs are defined as the AGN that fall in the UNCOVER color selection (shaded red area). In the left panels we report the fraction of cLRDs over AGN. The expected fraction is between 10 and 20%. In the right panels we report the fraction of cLRDs with MBH >  3 × 10−2Mstar) and with fEdd >  0.3. The yellow points report the colors considering only the AGN emission.

Keeping in mind the limitations of the simulation and of our model we have just discussed, the fact that the UNCOVER color selection, when applied to our model, is not picking galaxies, means that the UNCOVER color selection is well-suited to identify AGN, although it does not select the full AGN population: it recovers between 5% and 26% of all AGN with bolometric luminosity L_bol >  10⁴⁴ erg s⁻¹, depending on the assumed relation of MBH to stellar mass, AGN SED and dust properties, as shown in Fig. 9. Based on Harikane et al. (2023), the red population accounts for ∼20% of the broad-line selected AGN, but we note that differing color selections were involved in spectroscopic targeting for different surveys. We discuss this further in Sect. 5, but qualitatively, the selection misses AGN powered by MBHs with M_BH-M_star ratio < 3 × 10⁻², low f_Edd, as well as AGN with very high or very low column density. For the SMC-like dust case, a way to increase the completeness of the AGN selection is to extend the selection further red in m115W − m200W to −0.5 <  m115W − m200W <  3.0, leaving the rest of the selection criteria unchanged. Depending on the specific SED and dust model, the fraction of AGN in the selection increases by 50–100%. About half of the systems with 1 < m115W-m200W < 3 (plus the other UNCOVER cuts) are cases where the AGN already contributes at F200W, because of low column density (for the AGN) and/or low SFR (for the galaxy) and/or very high M_BH and f_Edd that lead the AGN to dominate everywhere. The others are galaxies with steep spectrum and often low SFR. MBHs with low M_BH-M_star ratio and/or f_Edd typically fail to enter the selection because they are unable to contribute sufficiently to the emission, which is dominated by the galaxy starlight. This was already noted in Volonteri et al. (2017, 2023) for models that did not include dust. In the “overm” case, the intrinsic MBH population is already skewed toward large MBH-to-galaxy ratios, and therefore it is easier for MBHs to contribute significantly to, or dominate, the emerging flux, as can be appreciated in Fig. 10 (right). In the “sEdd” case, instead, the redder peak of the SED increases the AGN flux in the red bands (redward of ∼F227W), easing their inclusion in a red selection (Fig. 10, left).

Fig. 10. Magnitudes for the “sEdd” and “overm” models and SMC-like dust. The orange triangles and gray circles show separately the flux of AGN and galaxies, when selecting systems with Lbol >  1044 erg s−1 and m444W < 28 (total flux, summing both AGN and galaxy contribution). The red dots have both of the above, plus the UNCOVER color selection, and show the combined galaxy+AGN flux. We have limited the y axis to 50 to highlight the relevant parameter space. At 115 and 200 the galaxy is almost always brighter than the AGN, whose contribution increases as the wavelength increases. For the “overm” case some AGN are brighter than the galaxy also at short wavelenghts.

The MW-like dust model³ generally leads to a smaller fraction of cLRDs than SMC-like dust, except in the “overm” case where the fraction instead increases. This is because fewer AGN have m115W − m200W >  1. Furthermore, with MW-like dust galaxies and especially AGN are fainter, therefore more are missed when selecting systems with m444W < 28. The number of AGN at the denominator decreases significantly, thus increasing the fraction of cLRDs over the full AGN population. Finally, with MW-like dust, the AGN with m115W − m200W >  1 are a minority, ∼3%.

Finally, we have assumed the same dust properties for the galaxy component and for the AGN. There is the possibility to have internal regions of a galaxy with a dust type (e.g., SMC) and external regions with another type (e.g., MW). Hence, the dust type could be of one kind for AGN and another kind for the galaxy component. If dust near AGN is SMC-like and dust in the rest of the galaxy is MW-like, cLRD selection is much more effective, and vice-versa. We may also expect variations in the extinction law for dust near AGN due to differential sublimation (smaller grains are more easily sublimated and have a larger sublimation radius). This generally translates into a flatter extinction curve (Gaskell et al. 2004; Gallerani et al. 2010; Antonucci 2023). We show an example in the Appendix, where we assume that dust for the stellar population is SMC-like and for the AGN is described by the empirical extinction curve from Gaskell et al. (2004).

5. Properties of AGNs and color-selected little red dots

We compare the properties of cLRDs to the full MBH population and AGN with L_bol >  10⁴⁴ erg s⁻¹ and m444W < 28 in Figs. 11 and 12. The clearest distinctive feature of cLRDs is their intermediate dust column densities. For high dust column densities the AGN is too buried and does not contribute sufficiently to the photometry. If the dust column density is too low, then AGN are not “red enough” as the dust-free case is not picked up by the selection, as shown in Fig. 8.

Fig. 11. Properties of color-selected cLRDs (red) with respect to all AGN with Lbol >  1044 erg s−1 (green) and with further m444W < 28 (blue) for SMC-like dust and model sEdd. A clear signature is the specific range in dust column density, Ndust, with respect to the global population. Further conditions are the elevated fEdd and ratio between MBH and Mstar.

Fig. 12. Properties of color-selected cLRDs (red) with respect to all AGN with Lbol >  1044 erg s−1 (green) and with further m444W < 28 (blue) for SMC-like dust and model overm. A clear signature is the specific range in dust column density, Ndust, with respect to the global population. Further conditions are the elevated fEdd and ratio between MBH and Mstar.

The bimodal behavior in the dust and hydrogen column density when imposing m444W < 28 (compare green and blue histograms) can be understood looking at Figs. 6 and 10. In galaxies with M_star < 10¹⁰ M_⊙, the intrinsic luminosities are smaller, therefore the magnitude cut selects galaxies (or lines of sight) with low column densities. In the most massive galaxies, the intrinsic luminosity is higher, therefore they can meet the magnitude cut even when column densities are higher: the high branch of the column density at M_star > 10¹⁰ M_⊙. This combines with the relation between MBH and galaxy luminosities, favoring overmassive MBHs with relatively high f_Edd. Even in the overm case, the median log(M_BH − M_star) ratio is higher for cLRDs than for all AGN with L_bol >  10⁴⁴ erg s⁻¹: −1.4 for cLRDs vis-a-vis −1.6 for all AGN with L_bol >  10⁴⁴ erg s⁻¹. For the sEdd case we have −2.4 for cLRDs, and −2.7 for L_bol >  10⁴⁴ erg s⁻¹. The median f_Edd = 1.0 is the same for cLRDs and L_bol >  10⁴⁴ erg s⁻¹ AGN in the sEdd case, while it’s −0.3 for cLRDs and −0.5 for L_bol >  10⁴⁴ erg s⁻¹ in the overm case.

In Fig. 6 we overlay in red the cLRDs to highlights how these points are interrelated. It is obviously easier for more massive galaxies to host more massive BHs, which at fixed f_Edd will more easily meet the bolometric luminosity requirement. However, if low-mass galaxies host sufficiently massive MBHs to power a bright AGN, they are “easier” hosts for cLRDs, because low-mass galaxies have lower column densities compared to their high-mass counterparts. The AGN can then more easily contribute to the total emission and skew colors. We also recall that in the “sEdd” case the redder intrinsic AGN SED helps to meet the color selection.

In terms of dust-to-gas ratio and SFR, color-selected AGN are evenly distributed within the global population. OBELISK galaxies are generally actively star-forming and have dust-to-gas ratios between 10⁻³ and 10⁻², and cLRDs span the same range, with a preference for low SFR. Certainly OBELISK galaxies have formed enough dust to explain the observed colors.

A question that often arises is what is the attenuation of cLRDs in terms of A_V. This is shown in Fig. 13 for the “overm” and “sEdd” cases. When considering only cLRDs, it is immediately obvious that they preferentially occupy a narrow range of A_V=1–5. This is because some attenuation is needed to meet the color selection, but enough AGN light has to escape in order to affect the colors. The “overm” case, which has the brightest AGN at fixed galaxy mass, allows for somewhat larger A_V(AGN), followed by the “sEdd” case. The median AGN A_V is 2.4, while for the “sEdd” case it’s 2.1. These values are somewhat higher than the median 1.6 found by Kokorev et al. (2024) and somewhat lower than the typical value of 3 found by Akins et al. (2024).

Fig. 13. AV for AGN versus AV stellar population, for AGN with Lbol >  1044 erg s−1 and m444W < 28 (first and third rows from top, where we have limited the y axis to AV = 50) and for systems selected as cLRDs (second and fourth rows from top). We only show SMC-like dust for the sEdd and overm models. AV for the stellar population of systems selected as cLRDs is moderately small with respect to the full extent of the distribution, while AV for the AGN is significantly skewed toward low values with respect to the global population.

Importantly, cLRDs are not the most attenuated sources in our analysis. Attenuation of the AGN component when we do not apply a magnitude cut can reach extreme values of 50 mag. These are cases that obviously preclude the AGN from showing any photometric signature. This is in contrast with the attenuation for the galaxy component, for which A_V is always limited to 3 mag. The median A_V for the stellar population is between 0.8 (sEdd) and 0.9 (overm). This is a consequence of the different (dust) column densities for the MBH and the stellar population. Fig. 6 highlights how NH for the stellar population, modulated by the location of young and old stars, is generally an order of magnitude smaller than for the centrally located MBH. The optical depth is proportional to the column density, and A_V is proportional to the optical depth. An order of magnitude larger column density therefore translates into an order of magnitude larger A_V.

In Fig. 14 (left) we turn to the global AGN population showing the AGN fraction as a function of stellar mass. The intrinsic active fraction, meaning the probability that a galaxy hosts an active MBH is unity in the model, but here we adopt an observer’s point of view, and consider the fraction of AGN bright enough to be detected, assuming that L_bol > 10⁴⁴ erg s⁻¹ is necessary to power sufficiently broad and bright emission lines. In all cases the fraction of AGN with L_bol > 10⁴⁴ erg s⁻¹ is an increasing function of M_star. This is because we have assumed a mass-independent f_Edd distribution, while MBH mass increases with galaxy mass. The “overm” case has significantly higher MBH masses, therefore the AGN fraction remains high also in low-mass galaxies. Adding the condition m444W < 28, set on the combined AGN+galaxy, decreases the AGN fraction in low-mass galaxies. This is because in galaxies with mass ≳3 × 10¹⁰ M_⊙ the stellar population by itself can have m444W < 28. This can be appreciated in Fig. 10: at the high-mass end cases with m444W of the AGN as faint at 40-50 can still meet the condition m444W < 28 on the total flux. Finally, the fraction of cLRDs is more sensitive to the MBH population properties (sEdd vs. overm; we do not show the distr case because too noisy) in galaxies with mass ≲10⁹ M_⊙. In the overm case it is easier for the brighter AGN to dominate on the galaxy colors.

Fig. 14. AGN populations. Left: AGN fraction, defined in terms of AGN bolometric luminosity as a function of galaxy mass. We recall that we have assumed an intrinsic active fraction of unity. Solid lines only require Lbol >  1044 erg s−1, dotted lines add m444W < 28 (galaxy and AGN combined) and dashed lines are for cLRDs (including the same bolometric and m444W cuts). Right: fraction of AGN-dominated sources at F444W as a function of m277W-m444W. Above m277W-m444W > 1.5 all sources are AGN-dominated.

In the right panel of Fig. 14 we touch on the observed compactness of sources. If we assumed that compactness meant AGN-dominated at F444W, then we find that sources with m277W-m444W > 1.5, have the AGN dominating photometry at F444W, for instance, m444W(AGN) < m444W(GAL). This agrees with the very high broad-line fractions and very high compact fractions found in observations for the same color cut. In other words, the redder the 277-444 color, the more AGN dominated is the light at F444W. We stress, however, that since we have no galaxy that without an AGN component has m277W-m444W > 1.5, and galaxies are found with such colors (Pérez-González et al. 2024), we cannot conclusively say that all sources with m277W-m444W > 1.5 are only AGN (see the discussion in Akins et al. 2024).

6. Dependence on the line of sight

The importance of the dust column density in determining whether an AGN is a cLRD leads us to explore the effect of observing a galaxy at fixed MBH and macro galaxy properties, through different lines of sight. In this section instead of drawing different M_BH and f_Edd for each line of sight, we use the same M_BH and f_Edd for all 12 lines of sight of a given galaxy.

In Fig. 15 we show examples of the photometry for a galaxy, where for different M_BH and f_Edd we show three different lines of sight. Fixing M_BH and f_Edd we can have, for the same source, an AGN-dominated case, a galaxy-dominated case and an intermediate case. The AGN dust column densities differ by 0.3 and 0.84 dex between the three lines of sight, while the effective dust column density for the stellar populations (Eq. (9)) differ at most by 0.1 dex. The line-of-sight difference is therefore much more pronounced for the AGN. This is because it is more sensitive to the position of few cells.

Fig. 15. Examples of photometry at z = 6 showing how different assumptions affect the emergent SEDs. For the same galaxy we show different MBH, fEdd and lines of sight. The galaxy photometry is calculated along the same line of sight as that of the AGN. The SFR is in units of M⊙ yr−1 and masses in units of M⊙ inside the logarithm. The NIRCAM response for the filters analyzed in this paper is shown in gray.

We have repeated the analysis shown in Fig. 9 and find the following cLRD fractions N(cLRDs)/N(AGN): 9% (distr), 13% (sEdd) and 18% (overm) for SMC-like dust and 4% (distr), 13% (sEdd) and 22% (overm) for MW-like dust, compatible with those obtained varying M_BH and f_Edd for each line of sight. As was discussed in the previous section the dust column density is as important as, if not more than, the intrinsic galaxy and MBH properties in determining whether an AGN can be classified as a cLRD.

It is instructive to see how at fixed galaxy and MBH properties different lines of sight move sources in the UNCOVER color selection. This is shown in Fig. 16 for 9 galaxies: three with mass < 10^8.5 M_⊙, four with mass 10^8.5 − 10^9.5 M_⊙ and two with mass > 10¹⁰ M_⊙. The results can be understood in light of Figs. 3, 15 and 8. When the AGN optical depth is high (large symbols) the AGN is too obscured and the source has colors typical of the stellar population, failing therefore to enter the UNCOVER selection being too blue in m277W-m444W and in m277W-m356W. When the AGN optical depth is low (small symbols) the AGN is unabsorbed, resulting also in colors that are too blue in m277W-m444W and in m277W-m356W. Only for intermediate optical depth (medium-size symbols) the AGN contribution blueward of F356W is suppressed resulting in redder colors and meeting the selection criteria.

Fig. 16. Colors for nine galaxies in which we explore how the appearance changes over 12 lines of sight, fixing MBH and fEdd. We show the sEdd and overm cases, for SMC-like dust. Each galaxy is shown with a different combination of color and symbol, therefore multiple identical markers indicate different lines of sight for the same galaxy+AGN intrinsic properties. The size of the symbol scales with the AGN optical depth. We show only cases with m444W < 28 (combined galaxy+AGN), but we do not impose any threshold on the AGN luminosity.

7. Conclusions

We have developed a comprehensive model of AGN photometry, including spatially resolved attenuation, in galaxies in the OBELISK simulation at z = 6. We first validate our model against observations, namely the AGN bolometric luminosity function and the X-ray obscured fraction. We then select AGN and cLRDs and investigate their properties. Our main results are as follows.

• Using realistic galaxy and AGN SEDs, with SMC- or MW-like dust we find that, as long as the AGN is sufficiently luminous with respect to the host galaxy, the AGN dominates in red NIRCAM filters and the galaxy in blue NIRCAM filters.

• AGN span a variety of gas and dust column densities, which depend both on the physical properties of the host galaxy and on the line of sight as fixed properties.

• The ISM contributes significantly to the total column density, and the contribution increases with galaxy mass and decreases with galaxy size.

• The simulation does not include any galaxy whose stellar population puts it in the little red dot UNCOVER color selection.

• Colors effectively select AGN with intermediate dust column density/attenuation. When attenuation is too low or too high AGN do not meet the m277W-m444W criterion, being too blue, and are not selected as cLRDs.

• cLRDs have high Eddington ratios and/or M_BH-M_star ratios. The host galaxies are within the general population, with only slightly lower SFR and not extreme dust-to-gas ratios.

• Modulating the SED at super-Eddington rates to mimic radiation trapping, thus making the SED redder, eases the selection of AGN as cLRDs. Alternatively, this means that cLRD selection favors super-Eddington accretors.

• cLRD selection is most sensitive to the MBH mass and accretion properties in galaxies with M_star <  10⁹ M_⊙.

• Sources with m277W-m444W > 1.5 are AGN-dominated, which can lead to such sources appearing compact, in agreement with observations.

• Given that only intermediate attenuation gives rise to cLRDs, there is a significant dependence on the line-of-sight properties, meaning that theoretical models require spatially-resolved galaxy properties.

Our analysis allows us to situate AGN selected as cLRDs within the AGN population, and understand their properties. Their spectral properties can be explained by AGN dominating in red filters and the host galaxy in blue filters because of intermediate attenuation and AGN that are intrinsically of luminosity comparable to that of the host galaxy.

Super-Eddington accretion, via trapping of UV radiation close to the inner edge of the accretion disk, can explain better the optical spectral shape of cLRDs (more power in red vs. blue) and there are arguments why it could be (but not always!) X-ray weak as many of the high-z AGN seem to be. In principle most of the high-z AGN could be super-Eddington. For cLRDs, super-Eddington helps pushing them toward the red selection, and therefore super-Eddington sources would naturally be more numerous with this selection. It is however unclear why purely spectroscopically selected AGN may also be predominantly super-Eddington. A possible explanation may be that super-Eddington accretion increases mass and luminosity of MBHs otherwise hosted in low-mass galaxies. Since low-mass galaxies are more numerous it would make it easier to preferentially select such sources (Trinca et al. 2022).

We have not found in the OBELISK simulation any galaxy that alone, without a AGN component, enters the cLRD selections. This is partly caused by the lack of emission lines in our model, which can boost photometry in the filter where they fall.

Based on these findings, future models need to account for spatially resolved line-of-sight attenuation properties, and model AGN SEDs with the same attention paid to physical properties as is currently done for galaxies. It is also important to understand if and how galaxies can enter the cLRD selection, and how common such galaxies are in theoretical models.

Data availability

All figures related to MW-like dust can be found at https://doi.org/10.5281/zenodo.14628499.

Acknowledgments

MV thanks Christina Eilers for clarifying the definition and usage of “little red dots”. MT and PD acknowledges support from the NWO grant 0.16.VIDI.189.162 (“ODIN”) and PD acknowledges support from European Commission’s and University of Groningen’s CO-FUND Rosalind Franklin program. RS acknowledges support from the PRIN 2022 MUR project 2022CB3PJ3 – First Light And Galaxy aSsembly (FLAGS) funded by the European Union – Next Generation EU. This research was supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP). This work has made extensive use of the Infinity Cluster hosted by Institut d’Astrophysique de Paris; we warmly thank Stéphane Rouberol for running smoothly this cluster for us. We acknowledge PRACE for awarding us access to Joliot Curie at GENCI@CEA, France, which was used to run most of the simulations presented in this work. Additionally, this work was granted access to the HPC resources of CINES under allocations A0040406955 and A0040407637 made by GENCI.

References

Appendix A: Additional cases

In this Appendix we report how results differ when (i) using different a color selection, or (ii) using an empirical dust law derived for AGN. We considered how the results differ when considering an alternative color selection of cLRDs following, for instance, Pérez-González et al. (2024), who adopt m277W−m444W> 1 and m150W−m200W< 0.5, plus m444W< 28. We show the results of our models for this selection in Fig A.1. This selection is similarly effective, especially in the SMC-like dust case, where more AGN have m150W−m200W around zero.

In Fig. A.2 we show results using the UNCOVER selection, where we assume that dust for the stellar population is SMC-like and for the AGN is instead described by the empirical extinction curve from (Gaskell et al. 2004). Whereas Li et al. (2024) find that this extinction curve provides a better fit for little red dots, we find instead that using this curve decreases the fraction of cLRDs among AGN (see also Ma et al. 2024b, for a detailed analysis on the spectrum of an AGN/cLRD). This is because Li et al. (2024) consider only the AGN contribution, while we include the galaxy contribution as well. In our case it’s the stellar population that provides the blue excess, while in the AGN-only case the flatter extinction curve of Gaskell et al. (2004) allows one to obtain a blue excess from the AGN only. In fact, if we remove the galaxy contribution to the photometry we find that the fraction of cLRDs among AGN increases to about 10-15% when using this extinction curve (compared to 4-8% including the galaxy component, as indicated in the figure).

Fig. A.1. SMILES color selection (shaded red area), with further Lbol >  1044 erg s−1 and m444W < 28 (combined galaxy+AGN). The AGN population is defined as meeting only the Lbol and m444W criteria. The yellow points report the colors considering only the AGN emission. The gray crosses highlight the sources that meet the UNCOVER cLRD selection.

Fig. A.2. UNCOVER color selection (shaded red area), with further Lbol >  1044 erg s−1 and m444W < 28 (combined galaxy+AGN). The AGN population is defined as meeting only the Lbol and m444W criteria. The yellow points report the colors considering only the AGN emission. Here we use the SMC-like extinction curve for the stellar population and the empirical extinction law derived in radio-loud AGN (Gaskell et al. 2004) for the AGN.

All Tables

All Figures

Fig. 1. Properties of simulated galaxies with Mstar >  108 M⊙ at z = 6, including galaxies that do not host MBHs. From top to bottom: star formation rate averaged over 10 Myrs in black, and over 100 Myrs in gray; mass-weighted stellar age, and mass-weighted stellar metallicity, where solar metallicity is assumed to be 0.02, a dust mass within twice the effective radius. Reff is defined as the half-mass projected radius in the plane orthogonal to the stellar angular momentum, computed as the (2D) radius such that 50% of the mass is within that radius. The violet line in the top panel shows the fit from Speagle et al. (2014).

In the text

Fig. 2. Sketch of our approach to modeling gas and dust absorption, showing how differential attenuation affects different parts of a given galaxy and lines of sight. We highlight two lines of sight with rays shown in purple and red arrows for each emitter. Darker gray colors show gas and dust clouds with larger optical depth. Stars’ colors indicate age: for instance young (blue) stars are generally still embedded in their birth clouds, while this is not necessarily the case for older (yellow and orange) stars. Since at least some star-forming regions are dusty, their contribution to the SED will be more affected by the extinction curve than other parts of the galaxy. The position of the AGN is marked in black, with thicker arrows. The example distribution of emitters and absorbers shown here has been chosen to highlight how the same source can be AGN- or stellar-dominated, depending on the line of sight. In this example, along the cumulative purple line of sight the AGN dominates over starlight: the AGN is seen through a clear line of sight, while much of the emission from the stellar population is absorbed. Vice-versa, along the cumulative red line of sight the AGN emission is sub-dominant to starlight.

In the text

	Fig. 3. Examples showing how different assumptions affect the emergent SEDs at z = 6. For the same galaxy we show different MBH, fEdd and NH, as well as different lines-of-sight for the emission from the stellar population. The NIRCAM response for the filters analyzed in this paper is shown in gray.
In the text

Fig. 4. Bolometric luminosity function for the actual MBHs in the simulation (“true”) and for the models described in Table 1, where “resc” means we have corrected the number densities for the overdensity bias: OBELISK is a protocluster. To obtain the cosmic average we have rescaled the OBELISK galaxy mass function to that of the NEWHORIZON simulation, which simulates an average region of the Universe. This figure assumes an active fraction of unity, and therefore represents an upper limit to the LF for each of the models. The results are compared to the luminosity functions derived in Matthee et al. (2024, M+24) and Greene et al. (2024, G+24) from JWST data and to fits to the evolution of the bolometric luminosity function derived in Shen et al. (2020, S+20, pre-JWST).

In the text

Fig. 5. Relation between hydrogen column density in the ISM and galaxy mass and half-mass radius at z = 6. Top panel: column density obtained integrating the column density from the position of the MBH to the host galaxy virial radius (“max” model). Bottom panel: column density obtained excluding the inner 80 pc around the MBH (“outmax80” model). The orange points show the mean column density for each galaxy, with size scaling with Mstar. For a random subset of 24 galaxies we show the column density for 12 different lines of sight, with the grayscale proportional to Mstar, as shown in the color bar.

In the text

Fig. 6. Hydrogen column density as a function of galaxy stellar mass. Top three rows: distribution of “effective” hydrogen column density for starlight obtained from the ratio of intrinsic and attenuated fluxes in each JWST filter we have considered (Eq. (9)). Bottom row: distribution of hydrogen column density seen by the AGN as a function of galaxy stellar mass at z = 6. This is the column density integrated above the dust sublimation radius, for example the column density that enters into calculating dust absorption. The red points highlight the cases that host a cLRD meeting the UNCOVER selection.

In the text

Fig. 7. Absorbed fraction, defined as the fraction of AGN with NH >  1023 cm−2 for various models compared to Vito et al. (2018) and Pouliasis et al. (2024). We include the 68% confidence region as shaded areas. The top panel shows the results including only the ISM gas, the bottom panel adds gas within the dust sublimation radius. For the cases including gas within the dust sublimation radius, we show the case with εr = 0.1, including spin evolution the absorbed fraction decreases to 60–70% (distr and sEdd), and 40% (overm).

In the text

	Fig. 8. Colors for galaxies only for SMC- (top) and MW-like (middle) dust extinction curves, and for galaxies and AGN without including dust (bottom). Galaxies at z = 6 do not enter the UNCOVER color selection highlighted by the shaded red area, nor do AGN without inclusion of dust. We consider here the “overm” case, which leads to the largest AGN contribution to colors.
In the text

Fig. 9. Colors for SMC-like dust. AGN are defined as sources having Lbol >  1044 erg s−1 and m444W < 28 (combined galaxy+AGN). cLRDs are defined as the AGN that fall in the UNCOVER color selection (shaded red area). In the left panels we report the fraction of cLRDs over AGN. The expected fraction is between 10 and 20%. In the right panels we report the fraction of cLRDs with MBH >  3 × 10−2Mstar) and with fEdd >  0.3. The yellow points report the colors considering only the AGN emission.

In the text

Fig. 10. Magnitudes for the “sEdd” and “overm” models and SMC-like dust. The orange triangles and gray circles show separately the flux of AGN and galaxies, when selecting systems with Lbol >  1044 erg s−1 and m444W < 28 (total flux, summing both AGN and galaxy contribution). The red dots have both of the above, plus the UNCOVER color selection, and show the combined galaxy+AGN flux. We have limited the y axis to 50 to highlight the relevant parameter space. At 115 and 200 the galaxy is almost always brighter than the AGN, whose contribution increases as the wavelength increases. For the “overm” case some AGN are brighter than the galaxy also at short wavelenghts.

In the text

	Fig. 11. Properties of color-selected cLRDs (red) with respect to all AGN with Lbol >  1044 erg s−1 (green) and with further m444W < 28 (blue) for SMC-like dust and model sEdd. A clear signature is the specific range in dust column density, Ndust, with respect to the global population. Further conditions are the elevated fEdd and ratio between MBH and Mstar.
In the text

	Fig. 12. Properties of color-selected cLRDs (red) with respect to all AGN with Lbol >  1044 erg s−1 (green) and with further m444W < 28 (blue) for SMC-like dust and model overm. A clear signature is the specific range in dust column density, Ndust, with respect to the global population. Further conditions are the elevated fEdd and ratio between MBH and Mstar.
In the text

Fig. 13. AV for AGN versus AV stellar population, for AGN with Lbol >  1044 erg s−1 and m444W < 28 (first and third rows from top, where we have limited the y axis to AV = 50) and for systems selected as cLRDs (second and fourth rows from top). We only show SMC-like dust for the sEdd and overm models. AV for the stellar population of systems selected as cLRDs is moderately small with respect to the full extent of the distribution, while AV for the AGN is significantly skewed toward low values with respect to the global population.

In the text

Fig. 14. AGN populations. Left: AGN fraction, defined in terms of AGN bolometric luminosity as a function of galaxy mass. We recall that we have assumed an intrinsic active fraction of unity. Solid lines only require Lbol >  1044 erg s−1, dotted lines add m444W < 28 (galaxy and AGN combined) and dashed lines are for cLRDs (including the same bolometric and m444W cuts). Right: fraction of AGN-dominated sources at F444W as a function of m277W-m444W. Above m277W-m444W > 1.5 all sources are AGN-dominated.

In the text

Fig. 15. Examples of photometry at z = 6 showing how different assumptions affect the emergent SEDs. For the same galaxy we show different MBH, fEdd and lines of sight. The galaxy photometry is calculated along the same line of sight as that of the AGN. The SFR is in units of M⊙ yr−1 and masses in units of M⊙ inside the logarithm. The NIRCAM response for the filters analyzed in this paper is shown in gray.

In the text

Fig. 16. Colors for nine galaxies in which we explore how the appearance changes over 12 lines of sight, fixing MBH and fEdd. We show the sEdd and overm cases, for SMC-like dust. Each galaxy is shown with a different combination of color and symbol, therefore multiple identical markers indicate different lines of sight for the same galaxy+AGN intrinsic properties. The size of the symbol scales with the AGN optical depth. We show only cases with m444W < 28 (combined galaxy+AGN), but we do not impose any threshold on the AGN luminosity.

In the text

	Fig. A.1. SMILES color selection (shaded red area), with further Lbol >  1044 erg s−1 and m444W < 28 (combined galaxy+AGN). The AGN population is defined as meeting only the Lbol and m444W criteria. The yellow points report the colors considering only the AGN emission. The gray crosses highlight the sources that meet the UNCOVER cLRD selection.
In the text

Fig. A.2. UNCOVER color selection (shaded red area), with further Lbol >  1044 erg s−1 and m444W < 28 (combined galaxy+AGN). The AGN population is defined as meeting only the Lbol and m444W criteria. The yellow points report the colors considering only the AGN emission. Here we use the SMC-like extinction curve for the stellar population and the empirical extinction law derived in radio-loud AGN (Gaskell et al. 2004) for the AGN.

In the text