© The Authors 2024

1 Introduction

In the upcoming years, ‘full-sky’ space and ground-based surveys, such as Euclid (Euclid Collaboration 2024g), MOONS (Cabral et al. 2020), Rubin/LSST (Ivezić et al. 2019), or eROSITA (Merloni et al. 2024), among others, will survey unprecedentedly large areas of the sky, gathering photometric and spectroscopic data for billions of galaxies and (at least) millions of active galactic nuclei (AGN). Synthetic data reproducing observed properties of astrophysical sources have become essential to enhance the scientific return from this new data (see e.g. Georgakakis et al. 2019, 2020; Comparat et al. 2020; Allevato et al. 2021; Bisigello et al. 2021; Euclid Collaboration 2024h). They are needed before the start of observations; for example, to define the selection bias, test data analysis pipelines, and develop and calibrate models and parameters for future analyses. Synthetic data can also help in understanding the observations; in particular, it can be used to estimate incompleteness and biases and to refine and validate hypotheses on the basis of the simulations.

Different methods can be used to generate mock catalogues (see Wechsler & Tinker 2018, for a review), ranging from physics-oriented models (hydrodynamic simulations or semi-analytic models) to data-oriented simulations (empirical models adopting observed scaling relations). Many of these simulations are primarily focussed on cosmological purposes, such as performing idealised tests on key cosmological probes like clustering and weak lensing (e.g. Korytov et al. 2019; Ishiyama et al. 2021; Hadzhiyska et al. 2023). However, to generate reliable forecasts from these measurements, it is crucial to incorporate realistic galaxy and AGN properties to accurately account for observational selection and systematic effects.

Recently, large-scale simulations that accurately reproduce galaxy properties have been developed (e.g. Kovacs et al. 2022; Wechsler et al. 2022; Gu et al. 2024; Euclid Collaboration 2024b), often adopting semi-empirical or hybrid methodologies. These methods not only save computation time but also ensure that results align well with most observed scaling relations. In addition, thanks to their optimised computational efficiency, they can also be quickly adapted to implement new empirical relations following new discoveries. Our work is conducted within this framework: we have adopted a semi-empirical method to produce realistic samples of synthetic galaxies and AGN, a component often neglected in previous similar works, necessary to enhance the realism of our simulation. At the same time, we preserve the link with the cosmic web traced by dark matter (DM) haloes: this connection is fundamental for deriving cosmological forecasts and linking the visible properties of galaxies and AGN to the DM distribution.

The inclusion of AGN in these catalogues is of the utmost importance; supermassive black holes (SMBHs) are believed to exist at the cores of most galaxies, and their masses are known to correlate closely with the masses, velocity dispersions, and luminosities of the host galaxy bulges (Richstone et al. 1998; Gebhardt et al. 2000; Merritt & Ferrarese 2001; Kormendy & Ho 2013). Furthermore, galaxies hosting actively accreting SMBHs – AGN – frequently present characteristic observational features (both in photometric and spectroscopic observations) across the full electromagnetic spectrum. Reproducing these features in realistic empirical catalogues is essential.

Although the aim of this work is to produce mock nearinfrared (NIR)/optical catalogues, we decided to implement AGN starting from X-ray-selected samples due to the efficiency of this selection in identifying AGN across a wide range of redshifts and luminosities, including those with high levels of nuclear obscuration (22 < log N_H/cm^–2 < 24) that may be overlooked in other bands (e.g. Padovani et al. 2017). Since X-ray AGN photons are generated through non-thermal processes near the central black hole, they serve as an ideal tracer of AGN activity (Mushotzky et al. 1993; Brandt & Hasinger 2005). Additionally, the X-ray band, especially for the brightest sources, experiences minimal contamination from the host galaxy, which predominantly emits at different wavelengths.

The structure of this paper is the following. In Sect. 2, we describe the methodology followed to create the galaxy mock catalogue that was used as a base for the consequent AGN mock. In Sect. 3, we describe in detail the methodology followed to populate the galaxy mock with AGN. Then, in Sect. 4, we validate the catalogue by comparing its outputs with observations that were not used for its calibration, and finally, in Sect. 5, we show one example of how this catalogue can be used to make predictions for future surveys, focussing on Euclid as a representative case.

Throughout this paper, a standard cosmology (Ω_m = 0.3, Ω_λ = 0.7 and H₀ = 70km s^–1 Mpc^–1) has been assumed. The stellar masses are given in units of solar masses for a Chabrier initial mass function (Chabrier 2003).

2 Mock galaxy catalogue: MAMBO

MAMBO (Mocks with Abundance Matching in BOlogna) is a workflow designed to construct an empirical mock catalogue of galaxies that can reproduce with accuracy their physical properties and observables, such as rest-frame and observed magnitudes and spectral features. A detailed description of the methodology and a validation of the galaxy properties can be found in Girelli (2021) and is briefly summarised in this section. A version of MAMBO, prior to the inclusion of AGN, has been used in Euclid Collaboration (2024d) to forecast the reliability of Euclid algorithms at recovering redshifts and physical parameters (stellar masses and star-formation rates) from Euclid photometry. In the rest of this paper, we focus our attention on the inclusion of AGN into this workflow. A schematic view of the steps explained in this section, as well as in Sect. 3, is presented in the flowchart of Fig. 1.

MAMBO takes in as inputs a few quantities from cosmological N-body DM simulations: the DM halo mass, expressed as M₂₀₀¹ for main haloes and Minfall² for subhaloes and orphans³; and the redshift, z, of each halo and subhalo. The sky coordinates RA and Dec are not necessary inputs, but can be used to reconstruct the cosmic web and derive clustering and environmental properties of different classes of objects.

The results presented in the following are based on a lightcone built by Henriques et al. (2015) using the Mock Map Facility (MoMaF, Blaizot et al. 2005) on the Millennium Simulation (Springel et al. 2005); namely, lightcone number 23, which was chosen because it presents a mass function that is the closest to the mean of all the 25 available lightcones. The lightcone spans from z = 0 to z = 10 and contains DM haloes with a minimum mass of M_halo ≳ 10^10.24 M_⊙ h^–1, which corresponds to 20 DM particles, and it covers an area of 3.14 deg². However, the method can be applied to any simulated catalogue of DM haloes and subhaloes.

In the first step, a galaxy with stellar mass ℳ is assigned to each DM halo by means of a stellar-to-halo mass relation (SHMR). The SHMR was derived using a subhalo abundance matching technique, and calibrated on the Millennium lightcones by means of observed stellar mass functions (SMFs) on the SDSS (York et al. 2000), COSMOS (McCracken et al. 2012) and CANDELS (Grogin et al. 2011) fields. A detailed description of the SHMR can be found in Girelli et al. (2020) and Girelli (2021).

Additionally, every galaxy in the mock is classified as pas- sive/quiescent (Q) or star-forming (SF) in a probabilistic way, following the relative ratio of the blue and red populations in observed SMFs. For this, we used the following SMFs: at z~0, the SMF evaluated by Peng et al. (2010) on the SDSS survey and divided into passive and SF using the rest-frame (U – B) colour; at 0.2 < z < 4, the SMF by Ilbert et al. (2013), derived on the COSMOS field and classified into red/blue using the rest frame colour selection (NUV – r) vs (r – J) (Ilbert et al. 2010).

At z ≥ 4, SMFs divided by SF/Q type are not available. Recent studies have tried to put constraints on this fraction (Merlin et al. 2017; Girelli et al. 2019; Xie et al. 2024), but still with large uncertainties. For this reason, the SF fraction was extrapolated from the results at lower redshifts, arriving at a maximum of f_SF = 99% at z = 6, which is kept constant at higher redshifts. This choice is motivated by the fact that at that redshift, the Universe is supposed to be too young (0.6 Gyr at z = 6) to contain any considerable fraction of quiescent galaxies. Recently, new results from JWST data are starting to find small samples of quiescent galaxies at 3 < z < 6 (Carnall et al. 2023, 2024; Alberts et al. 2023), while Looser et al. (2024) reported the finding of a quiescent galaxy at z = 7.3. We plan to revise these assumptions once sufficient data becomes available to accurately constrain the statistics of this population.

We show in Fig. 2 the stellar mass and redshift distributions resulting from applying the SHMR and the SF/Q classification process described above. As it is visible from this figure, the galaxy catalogue is complete at least down to ℳ~10^7.5 M_⊙, which is a consequence of the mass completeness of the DM halo catalogue. Additionally, it can be seen from this figure that the stellar mass distribution of quiescent galaxies follows the shape of a double Schechter function, as is often the case at redshifts where it is possible to observe also low mass galaxies (Drory et al. 2009; Peng et al. 2010; Ilbert et al. 2013; Davidzon et al. 2017; Weaver et al. 2023).

As a final step, other physical properties, as well as the photometry and spectra of the galaxies in the catalogue are retrieved with a modified version of the public code EGG (Empirical Galaxy Generator, Schreiber et al. 2017). EGG is a C++ code designed to generate an empirical mock catalogue of galaxies with realistic physical properties (such as the star formation rate, size, dust extinction, velocity dispersion, and emission line luminosities), where every galaxy is treated as a two component system (bulge + disc). Additionally, the code produces the observed and rest-frame photometry in any desired band, and the redshifted (observer-frame) spectra from the ultraviolet (UV) to the submillimeter. The code has been calibrated purely by using empirical relations to produce realistic observable properties. In order to assign an SED (spectral energy distribution) to each galaxy, EGG selects an optical template from a prebuilt library of Bruzual & Charlot (2003) models, covering uniformly the observed part of the UVJ plane (Williams et al. 2009), which separates quiescent from SF galaxies. Infrared SEDs are instead derived from a set of libraries (Chary & Elbaz 2001; Magdis et al. 2012; Schreiber et al. 2017) aimed at reproducing dust emission, characterised by the values of infrared luminosity of dust and PAHs (polycyclic aromatic hydrocarbons) at λ = [8, 1000] µm, the dust temperature, and the ratio of IR to 8 µm luminosity (IR8, Elbaz et al. 2011) for dust and PAH.

Afterwards, Gaussian emission lines with a given velocity dispersion⁴ were assigned to the SEDs of the bulge and the disc components. EGG can take as an input the stellar mass, ℳ, redshift, z, and type (SF or quiescent) of each galaxy, or produce these quantities by randomly extracting z and ℳ from the observed galaxy SMFs. We used EGG in the former configuration. More details about the way EGG works and about the modifications we made to the code are given in Sect. 3.5. It is worth stressing that the full MAMBO workflow can be applied to any DM simulation, and, similarly, the method presented in Sect. 3 to populate galaxies with AGN can be applied to any mock galaxy catalogue containing information about galaxy stellar mass, redshift, and galaxy type.

Fig. 1 Schematic view of the full workflow presented in this paper to create the galaxy and AGN mock catalogue. Yellow boxes represent the result of each step, while blue boxes show the necessary inputs. Steps within the dotted green box take place within the modified version of the public code EGG. A detailed description is given in Sects. 2 and 3.

Fig. 2 Distribution of stellar mass and redshift of all the galaxies in the MAMBO lightcone used in this work, separated into the SF (blue) and quiescent (red) populations. The dashed black line shows the total (SF+Q) galaxy population.

3 Adding AGN to the catalogue

In this section, we present the methodology adopted to populate our galaxy catalogue with AGN, which uses a completely empirical methodology comprising the following steps: i) we flagged each galaxy from the lightcone as hosting an AGN or not, following a probabilistic method that depends on ℳ and z (Sect. 3.2); ii) every object flagged as AGN was assigned an intrinsic X-ray luminosity (Sect. 3.3); iii) we separated the AGN population into optically unobscured or obscured AGN (type 1 or type 2, respectively) by means of their intrinsic X-ray luminosity (Sect. 3.4); iv) we built the observed spectra and photometry of both type 1 and type 2 AGN with the help of photoionisation models of AGN narrow-line regions and parametric SEDs of typical quasi-stellar objects (QSOs), which model both the continuum and emission lines (see Sect. 3.5). A schematic view of this workflow is given in Fig. 1.

3.1 Specific accretion rate

The specific accretion rate (λ_SAR ∝ L_X/ℳ)⁵ is a directly measurable quantity that has been extensively used in the literature (Bongiorno et al. 2012; Aird et al. 2012; Georgakakis et al. 2014) as a proxy for the rate of accretion onto the SMBH (Ṁ_BH) relative to the stellar mass of the host galaxy. Furthermore, if a bolometric correction and an M_BH/ℳ scaling relation are assumed, the specific accretion rate can also be regarded as a proxy for the Eddington ratio of the SMBH, λ_Edd = L_bol/M_BH (Bongiorno et al. 2016; Aird et al. 2018).

In the following subsections, we use two specific accretion rate distribution functions (SARDFs) with different definitions. On the one hand, Bongiorno et al. (2016) defines the specific accretion rate as λ_SAR = L_X/ℳ. On the other hand, Aird et al. (2018) defines the specific black hole accretion rate (λ_sBHAR) as the dimensionless quantity: $${\lambda }_{\text{sBHAR}}=\frac{k_{\text{bol}}{L}_{\text{X}}}{1.3\times {10}^{38}\text{erg} {\text{s}}^{-1}\times 0.002\frac{ℳ}{M_{\odot }}},$$(1)

where k_bol is a bolometric correction factor (L_bol = k_bol L_X), assumed to have a constant value of k_bol = 25. This definition also assumes a constant scaling relation of M_BH = 0.002ℳ (Marconi & Hunt 2003, assuming also ℳ ≈ M_bul). Under these assumptions, λ_sBHAR = λ_Edd, and therefore an AGN accreting at 1% of the Eddington limit, would have λ_sBHAR~10^–2. It is possible to convert between the two definitions of the specific accretion rate as λ_SAR ≈ 10³⁴ λ_sBHAR. Therefore, with the same assumptions, the 1% of the Eddington limit corresponds to ${\lambda }_{\text{SAR}}={10}^{32}{\text{erg s}}^{-1}{M}_{\odot }^{-1}.$ We note that both in Bongiorno et al. (2016) and in Aird et al. (2018), this was set as the lower limit to define a galaxy as hosting an AGN.

3.2 Fraction of AGN

In order to derive the probability of a galaxy with a given stellar mass, ℳ, and at a given redshift, z, hosting an active nucleus, p (AGN | ℳ, z), we made use of the AGN host galaxy mass function (HGMF) derived by Bongiorno et al. (2016, hereafter B16). In B16, the authors studied a sample of 877 hard (2–10 keV) X- ray selected AGN from the XMM-COSMOS point-like source catalogue (Hasinger et al. 2007; Cappelluti et al. 2009; Brusa et al. 2010; Bongiorno et al. 2012) in the redshift range 0.3 < z < 2.5 and with a limiting flux of F_X~3 × 10^–15 erg s^–1 cm^–2. This sample was also selected with a stellar mass limit of ℳ > 10^9.5 M_⊙ and a specific accretion rate of ${\lambda }_{\text{SAR}}={10}^{32}{\text{erg s}}^{-1}{\text{M}}_{\odot }^{-1}.$

To derive the AGN HGMF and the SARDF, B16 corrected for the stellar mass incompleteness of the sample down to the above-mentioned limit in λ_SAR. Additionally, B16 also accounted for the incompleteness due to the sources that were missed in the sample because of their high column density, N_H. This was done considering column density values in the range of 20 < log N_H/cm^–2 < 24, and therefore not including Compton-thick AGN, which are heavily obscured objects with log N_H/cm^–2 > 24. Although a significant fraction of AGN are expected to be CTK sources, their exact fraction is still a matter of debate (e.g. Buchner et al. 2015 found a constant fraction of f_CTK∼35%, independent of redshift and accretion luminosity, while Pouliasis et al. 2024 found a much lower fraction of f_CTK ∼ 17% for 3 ≲ z ≲ 6 AGN). Furthermore, the redshift and luminositydependence of this fraction is yet not fully understood (e.g. Ricci et al. 2017a). Therefore, the inclusion of these sources is left to future work.

It is worth noting that the choice of a minimum value of λ_SAR sets the definition of AGN used in this work. In general, different criteria can be used to select AGN based on their X-ray emission. A common method is to select sources emitting above a threshold intrinsic X-ray luminosity (generally L_X > 10⁴² erg s^–1, e.g. Brandt & Alexander 2015). The AGN can also be selected according to the relative X-ray emission from non-nuclear mechanisms (such as X-ray binary stars and hot gas emission) with respect to the total L_X of the galaxy (e.g. Birchall et al. 2020). It is well known that selecting AGN with these different criteria can lead to significant differences in their observed fraction (see e.g. Birchall et al. 2022), which should be taken into consideration when using this catalogue.

In B16, the authors derived the HGMF by jointly fitting the SMF and the SARDF (see Sect. 3.3), with the X-ray luminosity function (XLF) as an additional constraint. For this purpose, they used a maximum likelihood method to determine the HGMF and the SARDF as a bivariate distribution function of the stellar mass and a specific accretion rate, Ψ(ℳ, λ_SAR, z). As a result of this approach, the HGMF cannot be expressed as a simple analytic function. Instead, the authors provide an analytic approximation by performing a least-squares fit to the HGMF with a Schechter (Schechter 1976) function: $$\text{Φ}(ℳ)\text{d}ℳ={\text{Φ}}^{\star }{\left(\frac{ℳ}{{ℳ}^{\star }}\right)}^{\alpha }\exp \left(-\frac{ℳ}{{ℳ}^{\star }}\right)\text{d}\left(\frac{ℳ}{{ℳ}^{\star }}\right),$$(2)

evaluated at the centre of three redshift bins (0.3 < z < 0.8, 0.8 < z < 1.5 and 1.5 < z < 2.5). The best-fit parameters of this fit are given in Table 2 of B16⁶.

The original slope (α) from the Schechter function in the first two redshift bins would produce an unrealistic overestimate of low-mass AGN when extrapolating this model to ℳ < 10^9.5 M_⊙, which is the stellar mass limit in the sample of B16, while our catalogue spans to lower stellar masses, as is visible in Fig. 2 (see Appendix A for an example using the original slopes). Therefore, we re-derived these quantities by fitting the 1 / V_max points shown in B16 with a Schechter function using Eq. (2), obtaining α = –0.25 ± 0.09 for 0.3 < z < 0.8, and α = –0.19 ± 0.11 for 0.8 < z < 1.5. For the highest redshift bin, our fit was compatible with the value shown in Table 2 of B16; therefore, we did not modify it.

We then calculated p (AGN | ℳ, z) as the ratio between the SMF of the galaxies of the MAMBO lightcone and the HGMF from B16. Because the MAMBO lightcone covers a wider range of redshifts (0 < z < 10), we interpolated and extrapolated p (AGN | ℳ, z) from the centre of each redshift bin. For this, we assumed a minimum fraction, p (AGN | z = 0) = 0.01, at all ℳ when extrapolating at z < 0.55, while we maintained a constant fraction when extrapolating at z > 2. Afterwards, every galaxy was statistically assigned as either a host of an AGN or not a host, with a Bernoulli trial proportional to p(AGN) (i.e. by comparing p(AGN) to a random number extracted from a uniform probability distribution from 0 to 1).

Although the choice of a minimal fraction of 1% at z ∼ 0 for all masses is a rough approximation, recent studies of the local Universe motivate this assumption. For example, Birchall et al. (2022) studied 917 X-ray selected AGN (found as XMM counterparts of 25 949 SDSS galaxies) with z ≤ 0.3 3, which corresponds to a global AGN fraction of 3.5%. Instead, when selecting AGN by means of their accretion rate they found a fraction of about 1%, constant over stellar masses of 8 < log(ℳ/M_⊙) < 12, and increasing from 1% to about 10% with redshift. In Williams et al. (2022) the authors studied 213 Chandra X-ray counterparts of 280 nearby (<120Mpc) galaxies from the Palomar sample and classified 14 (6.6%) of them as Seyferts, while only four (1.9%) of them have L_X > 10⁴². Similarly, Osorio-Clavijo et al. (2023) studied 138 Chandra X-ray counterparts of the CALIFA sample, with a wide range of stellar masses and z < 0.1 and found an AGN fraction of 5%. At z > 2, we chose to have a constant AGN fraction: despite being a rough approximation, with this choice we are able to reproduce the observed X-ray luminosity functions up to z = 7, as is shown in Sect. 4.2.

In the upper panels of Fig. 3, we show the best-fit Schechter function of the AGN HGMF (both the original from B16 and our modified fit) together with the SMF of the galaxies of the MAMBO lightcone and the AGN mass function of our catalogue, derived with the methodology described above. We also show the AGN HGMF computed using the V_max method, derived in B16, as an additional consistency check to the Schechter model. As was expected by construction, the AGN mass functions of MAMBO reproduce those of B16. An exception is the highest redshift bin, where the low-mass end of the HGMF of MAMBO is higher, but still compatible with the V_max estimate, due to the fact that we kept a constant p (AGN | ℳ, z) when extrapolating at z > 2 (see Fig. 4). In the lower panels of Fig. 3 we show the fraction of AGN over the galaxy population as a function of stellar mass and at a given redshift bin, p (AGN | ℳ, (z)). As was expected, this fraction increases with increasing stellar mass of the host galaxy.

We show in Fig. 4 the probability of a given galaxy hosting an AGN as a function of redshift and in different mass bins for the full lightcone, after performing the interpolations and extrapolations described above. Each dot in this figure corresponds to a different galaxy, which is flagged as an AGN (or not) in a statistical way depending on this probability. This probability, at all masses, increases with redshift, reaching a maximum around z ≃ 1 and decreasing at higher redshifts. For comparison, Fig. 4 also shows the fraction of AGN derived in Aird et al. (2018, hereafter A18), in which they derived the duty cycle of a sample of NIR-selected galaxies matched with X-ray data. We note that the fractions differ significantly, being higher in this work, in particular at ℳ < 10¹⁰ M_⊙. We show in Appendix A that using the AGN fractions from A18 as an input in our workflow produces an X-ray luminosity function that tends to be underestimated when comparing it to the observed ones, especially for z ≲ 1 and z ≳ 3.

Fig. 3 Stellar mass functions of galaxies and AGN in the three redshift bins studied in B16. The green and purple shaded areas correspond to the galaxies and AGN from the MAMBO lightcone, respectively. We also show the best-fit Schechter function to the HGMF of AGN from Bongiorno et al. (2016) with a dotted black line, and the modified Schechter fit used in this work with a dashed black line. The shaded grey area shows the uncertainty in a. Black dots with error bars show the AGN mass function computed using the Vmax method, derived in B16. The lower pannels show the AGN fraction derived as the ratio of the modified Schechter function and the galaxy mass function from MAMBO (dashed blue line), and the same but using the original Schechter function from B16 (dotted black line).

Fig. 4 Probability of a galaxy hosting an AGN as a function of redshift. Each point corresponds to a galaxy in the MAMBO lightcone, colour-coded in different ℳ bins indicated in the legend. For comparison, the points with error bars connected by dashed lines show the AGN fraction derived in A18, also in different mass bins. Vertical dotted lines mark the centre of the three redshift bins studied in B16; namely, z = 0.55, 1.15, 2.0.

3.3 X-ray luminosity

In order to assign to each object flagged as AGN an intrinsic luminosity in the hard ([2–10 keV]) X-ray band, we first assigned a specific accretion rate. At 0 < z < 2 we make use of the SARDF derived in B16. As is explained in Sect. 3.2, B16 derived the SARDF and the HGMF simultaneously as a bivariate distribution function of ℳ and λ_SAR; that is, Ψ(ℳ, λ_SAR, z). As was the case with the HGMF, this implies that the SARDF cannot be expressed as a simple analytic function of ℳ and z. Instead, we used an analytic approximation of the SARDF (evaluated at the centre of three redshift bins) described as a double power-law (DPL) of the form $$\text{Φ}\left({\lambda }_{\text{SAR}},ℳ\right)=\frac{{\text{Φ}}_{\lambda }^{\star }}{{\left(\frac{{\lambda }_{\text{SAR}}}{{\lambda }_{\text{SAR}}^{\star }(ℳ)}\right)}^{-{\gamma }_1}+{\left(\frac{{\lambda }_{\text{SAR}}}{{\lambda }_{\text{SAR}}^{\star }(ℳ)}\right)}^{-{\gamma }_2}},$$(3)

where the mass dependence is given by log ${\lambda }_{\text{SAR}}^{\star }(ℳ)=\log {\lambda }_{\text{SAR},0}^{\star }+k_{\lambda }\left(\log ℳ-\log {ℳ}_0\right)$, where $\log {\lambda }_{\text{SAR},0}^{\star }=33.8,\log {ℳ}_0=11$, and k_λ = 0.58. The best-fit values of the normalisation, ${\text{Φ}}_{\lambda }^{\star }$, and slopes, γ₁, γ₂, of the DPL evaluated at the centre of the three redshift bins are given in Table 3 of B16.

Using Eq. (3), we computed the SARDF at five values of ℳ, from log(ℳ/M_⊙) = 9.75 to log(ℳ/M_⊙) = 11.75 in steps of 0.5 dex in ℳ (Fig. 5). We then divided the AGN in our lightcone into five mass bins, each of them centred at one of the abovementioned ℳ and all of them with a width of 0.5 dex in ℳ, except for the lowest-mass bin, extending up to the lowest ℳ in our lightcone. Then, we normalised the λ_SAR distribution functions by dividing them by the integrated AGN density (Mpc^–3) in each bin of ℳ and z, and therefore converting these distributions into probability distributions. We assigned to each AGN a value of λ_SAR by randomly extracting values from the corresponding probability distribution at each ℳ and z bin.

At higher redshifts (z > 2), we instead made use of the accretion rate distributions from A18, who studied a sample of 1797 X-ray Chandra counterparts of 126 971 NIR-selected galaxies up to z ~ 4 and with stellar and 8.5 < log(ℳ/M_⊙) < 11.5. Using a methodology similar to that used for the B16 accretion rate distributions previously described, we constructed the cumulative distribution function relative to each p (log λ_sBHAR | ⊙, z) distributions in the range of –2 < log λ_sBHAR < 1 and 2 < z < 4, and used it to assign a value of λ_sBHAR to each object previously labelled as an AGN at z > 2. Since A18 calculated the λ_sBHAR distributions separately for SF and quiescent galaxies, we used both sets of distributions for the galaxies in our catalogue that are split into these two classes. We also explored the possibility of using the accretion rate distributions from A18 to assign L_X at all redshifts, but we decided to use the methodology presented in this section since it produces an XLF that is in better agreement with the observed ones (see Appendix A).

After assigning every AGN in the lightcone a specific accretion rate, we converted this into the intrinsic X-ray luminosity (in the 2–10 keV band). We show in Fig. 6 the histogram of the assigned LX in different redshift bins for the entire population of AGN. We note that, as is visible in Fig. 6, a large fraction of AGN have luminosities of L_X < 10⁴² erg s^–1; that is, below a classical threshold applied to define an X-ray emitter as an AGN. The reason for this is twofold: on the one hand, as was previously mentioned, we started from a sample of AGN selected above a given accretion rate limit and not above a L_X limit. On the other hand, the galaxy mock probes low stellar masses (ℳ ∼ 10⁷ M_⊙), which translates into low X-ray luminosities.

In Fig. 7, we show the distribution of the AGN in our mock in the L_X–ℳ plane. As was expected by construction, all the AGN have accretion rates above ${\lambda }_{\text{sBHAR}}>-0.01(\equiv {\lambda }_{\text{SAR}}>{10}^{32}\text{erg} {\text{s}}^{-1}{\text{M}}_{\odot }^{-1}$), and only a small fraction of them are above the Eddington limit of λ_sBHAR = 1.

Fig. 5 Specific accretion rate distribution function (SARDF) computed at different values of ℳ using Eq. (3), in the three redshift bins from B16.

Fig. 6 Distribution of intrinsic hard-band X-ray luminosity (LX) of all the AGN in the lightcone, separated into type 1 (blue) and type 2 (orange) AGN, in different redshift bins. We also show with dashed and solid lines the subpopulations of type 1 and 2 in SF and quiescent galaxies. The vertical dotted line marks LX = 1042 erg s–1, the threshold typically applied to separate X-ray emission from AGN or from other origins.

3.4 Obscuration model

The observed spectra of AGN in different bands can vary significantly depending on the level of obscuration of the nucleus by its surrounding material along the line of sight. In this context, AGN are generally separated into two main sub-classes, namely type 1 and type 2 AGN (unobscured and obscured, respectively). In the UV/optical/NIR bands, type 1 AGN are characterised by the presence of broad optical emission lines (full width at halfmaximum >1000 km s^-1) and a SED with a bluer continuum, while type 2 AGN typically do not present these characteristics. In reality, AGN display a much wider observational variety, which allows us to classify them in many other sub-classes (e.g. Antonucci 1993; Urry & Padovani 1995; Spinoglio & Fernández- Ontiveros 2021), but such classification is beyond the scope of this work.

In order to separate the AGN of our catalogue into optically obscured and unobscured, we make use of the results presented in Merloni et al. (2014). In this study, the authors studied a sample of 1310 AGN selected from the XMM-COSMOS point-like source catalogue (Hasinger et al. 2007; Cappelluti et al. 2009) with a limiting rest-frame X-ray flux of F_X = 2 × 10^–15 erg cm^–2 s^–1 in the redshift range 0.3 < z < 3.5. For this sample, the authors studied the luminosity dependence of the (optically) obscured fraction of AGN, and found a relation, almost redshift independent, of the form $$F_{\text{obs}}=A+\frac{1}{\pi }\mathrm{atan}\left(\frac{l_0-\log L_{\text{X}}}{{\sigma }_x}\right),$$(4)

where F_obs is the fraction of type 2 AGN with respect to the total population, L_X is the intrinsic X-ray luminosity in the hard band, and the best-fitting parameters are A = 0.56, l₀ = 43.89 and σ_x = 0.46. For this work, we assumed this relation to be valid for the full redshift range 0 < z < 10 without any redshift evolution. Using this relation, we assign to every AGN a probability of being optically obscured depending on their L_X, and statistically separate every AGN in type 1 and type 2.

Additionally, the X-ray emission from AGN (especially in the soft band) can be heavily obscured by gas and dust surrounding the SMBH. To quantify this effect, we assign to each AGN a value of absorption column density (N_H), following the absorption function presented in Ueda et al. (2014), within the range 20 < log N_H/cm^–2 < 24. The choice of this range is consistent with the fact that, as stated in Sect. 3.2, Compton-thick AGN (logN_H/cm^–2 > 24) are not included in our model. Following Sect. 3 of Ueda et al. (2014), we assign N_H in a probabilistic way as a function of L_X and redshift. This allows us to classify the AGN in our catalogue as X-ray type 1 and type 2, where type 1 X-ray AGN are defined as an object with log N_H/cm^–2 < 22, and vice-versa.

It is worth noting that in nature these two classifications, although they on average correlate with each other, do not match perfectly; that is, some objects present optical obscuration but not X-ray obscuration and vice versa (e.g. Marchesi et al. 2016b; Merloni et al. 2014). Because the aim of this work is to produce a catalogue to be used in optical/NIR surveys (e.g. Euclid, MOONS), in the following sections we focus on type 1 AGN defined as optically unobscured sources, regardless of their X-ray obscuration.

In Fig. 6, we show the LX distribution separated into type 1 and type 2 AGN. It is visible from this figure that most of the AGN with L_X < 10⁴² erg s^–1 are classified as type 2 AGN in our catalogue.

In Fig. 7, we show also the distribution of sources from the Chandra COSMOS Legacy Spectral Survey (Marchesi et al. 2016b) in the L_X–ℳ plane, selected above a minimum X-ray flux F_X = 1.9 × 10^–15erg s^–1 cm^–2, which is the flux limit reported in Marchesi et al. (2016b). We applied the same flux cut to the AGN from our catalogue (contours in Fig. 7) in order to compare with the COSMOS data. With this flux cut, our AGN and the COSMOS AGN are limited to z < 3.5. We observe a reasonable agreement between the MAMBO and the COSMOS distributions, although MAMBO type 1 AGN seem to cover a narrower range in stellar mass than the COSMOS ones, and MAMBO type 2 AGN have on average lower X-ray luminosities than their COSMOS counterparts.

In Fig. 8, we show the fraction of optically obscured AGN as a function of L_X for the presented lightcone, and compare it with the calibration used to derive it; that is, Eq. (4). We also show the fraction of objects that are obscured both in optical and X-ray, and those that are unobscured in both bands.

Fig. 7 AGN from our lightcone scattered in the LX–ℳ plane, separated into type 1 (blue) and type 2 (orange) AGN. The points in the central panel show the distribution of all the sources, while the solid contours and the histograms from the upper and right panel correspond to sources selected above FX = 1.9 × 10–15 erg s–1 cm–2. We show also the distribution of sources from the Chandra COSMOS Legacy Spectral Survey (Marchesi et al. 2016b), with the same cut in FX, and separated into type 1 (cyan) and type 2 (pink) AGN. The contour levels represen iso-density lines, corresponding to the 50th, 75th, 90th, and 99th percentiles of the distribution. The dash-dotted and dotted black lines mark the locus where λSBHAR = 0.01 and 1, respectively, which assuming a mean bolometric correction of kbol = 25 and a constant mass ratio of black hole to host galaxy of MBH ≈ 0.002ℳ correspond approximately to 1 and 100% of the Eddington limit, respectively.

Fig. 8 Fraction of optically obscured (type 2) AGN as a function of intrinsic X-ray luminosity. The solid green line shows the relation used to separate AGN (Eq. (4), originally from Merloni et al. 2014) into optically obscured/unobscured, while the dashed purple line shows the actual fraction of type 2 AGN in our catalogue (with the Poissonian uncertainty shown by the shaded area). The other two dashed lines give information about the X-ray obscuration (where type 1 X-ray AGN are defined as objects with logNH/cm–2 < 22 and vice-versa): the dashed orange line shows the fraction of objects classified as type 2 both in optical and the X-ray, while the dashed red line represents the fraction of all the AGN that are type 1 in both bands at the same time.

3.5 AGN emitted spectra

After every galaxy in the catalogue has been classified as either hosting an AGN or not, and AGN have been characterised in terms of their X-ray luminosity and optical obscuration, we employ the publicly available code EGG (Empirical Galaxy Generator, Schreiber et al. 2017) to assign the rest of physical properties and observables (e.g. bulge/disc ratio, dust attenuation, etc.). Additionally, EGG allows us to generate the photometry of every object in any desired band, as well as the complete SED from the UV to the submillimeter.

In EGG, each galaxy is represented as a two-component system composed of a disc and a bulge, each of which is associated with a distinct SED, selected from a predefined lookup table. The choice of the galaxy SED is based on specific recipes that are tied to three main galaxy properties: its total stellar mass, its redshift, and its type (SF or quiescent).

We modified the original code by adding a third component to account for the AGN emission. For type 1 AGN (optically unobscured), this component includes the continuum and narrow and broad line emission typical of QSO, while for type 2 it accounts only for narrow-line emission (i.e. the galaxy stellar continuum is assumed to dominate at all wavelengths). For normal galaxies, this component remains effectively null.

3.5.1 Type 1 spectral energy distributions

In order to construct the SEDs of type 1 AGN in the lightcone, we made use of the parametric model developed by Temple et al. (2021). This model allows for the generation of synthetic quasar SEDs over the rest-frame wavelength range 912 Å to 3 µm. These synthetic SEDs have been shown to accurately reproduce, to a high degree of accuracy, the observed-frame optical and NIR colours of large samples of quasars, over redshifts 0.2 ≤ z ≤ 7 and absolute magnitude in the i-band –29 < M_i < –22.

The observed variety in emission line properties is included in the model through the interpolation between two emissionline templates, which correspond to the observed limits of very strong and very weak emission in terms of the equivalent width of high-ionisation UV lines, such as CIV. Furthermore, observations from quasar spectra show that the equivalent width of strong emission lines is anti-correlated with the intrinsic luminosity of the source. This phenomenon is known as the Baldwin effect (Baldwin 1977). In order to reproduce this phenomenology, the model from Temple et al. (2021) incorporates a single parameter (emline_type) that allows for the generation of spectra with different emission line properties. This parameter is related to the absolute magnitude M_i of the quasar by means of the following empirical equation: $$\text{emline\_type\hspace{0.17em}}=0.183\times \left(M_i\text{+27}\right).$$(5)

For this work, we constructed the QSO SED library by changing only this parameter between its minimum and maximum values (–2 to +3), ranging from spectra with weak, highly blueshifted lines to those with strong, symmetric lines.

For each type 1 AGN, we modified the code EGG to select a SED from the pre-built library of QSO SEDs following Eq. (5). The SED is rescaled to the L_UV (at 2500 Å) that corresponds to that object according to the observed L_X–L_UV , as is presented in Bisogni et al. (2021). This relation is parametrised as $$\log \left(L_{\text{X}}\right)=\gamma \log \left(L_{\text{UV}}\right)+\beta ,$$(6)

where the proxies for the UV and X-ray emissions correspond to the monochromatic rest-frame 2500 Å and 2 keV luminosities, respectively. The authors in Bisogni et al. (2021) find an average value for the slope of the relation of γ = 0.58 ± 0.06 up to z ≈ 4.5 and a dispersion of δ = 0.24 dex, which we used also for this work. For this, we calculated the monochromatic luminosity at 2 keV from LX following the general relation between the total luminosity in the band (L_{E1 -E2}) and at a monochromatic energy (L_E); that is, $$L_E=\frac{(2-\text{Γ})E^{1-\text{Γ}}}{E_2^{2-\text{Γ}}-E_1^{2-\text{Γ}}}{L}_{E_1-E_2},$$(7)

where in this case E = E₁ = 2 keV, E₂ = 10 keV, and we used a fixed value for Γ = 1.8. Therefore, following Eqs. (6) and (7) we derive the monochromatic luminosity L₂₅₀₀ from the full band luminosity L_X, and use it to rescale the SED.

We apply the reddening by dust attenuation to the QSO SED using an empirically derived extinction curve presented in Temple et al. (2021). This extinction curve is similar to that of the Small Magellanic Cloud for λ ⪆ 1700 Å, while it increases less rapidly for shorter wavelengths. The E(B– V) of each source is chosen randomly from the distribution presented in Fig. 2 of Lusso et al. (2013), which peaks at low E(B – V), with a median value 〈E(B – V)〉 = 0.03.

Finally, the galaxy SED (bulge + disc components) is added to that of the QSO in order to create the composite spectrum. To quantify how much each type 1 AGN is dominated by the hostgalaxy stellar light or by the nuclear emission, we derived the quantity f_AGN , which is defined as $$f_{{\text{AGN}}_{\text{band}}}=\frac{F_{\text{band }}^{\text{AGN}}}{F_{\text{band }}},$$(8)

where $F_{\text{band }}^{\text{AGN}}$ is the flux in a given band from the AGN component only, and F_band is the equivalent for the total spectrum (i.e. disc + bulge + AGN). With this definition, f_AGN ranges from 0 to 1, where an object with f_AGN = 0(1) would be completely dominated by the galactic (AGN) component, and f_AGN = 0.5 corresponds to the limit where both the galactic and the nuclear components contribute equally to that specific band. We derived this quantity in two bands: using the rest-frame magnitude from the GALEX FUV filter (which corresponds to the wavelength range where AGN emission typically peaks), and using the observed magnitude m_H (see Table 1).

In Fig. 9, we show an example of rest-frame SED for both a type 1 and a type 2 AGN from the lightcone. We show also in the figure the value of the most relevant physical parameter related to the construction of the SED, such as M, LX , SFR and, only for the type 1 AGN, $f_{{\text{AGN}}_{\text{UV}}}$ and $f_{{\text{AGN}}_H}$.

Fig. 9 Some examples of SEDs of AGN from our lightcone. We show three representative cases: in the upper panel, we show a type 1 AGN with dominant QSO component in the rest-frame UV band, while in the middle panel, we show a type 1 AGN with similar QSO and host-galaxy contribution for the same band. The lower panel shows a type 2 AGN, with a zoom-in subpanel showing the Hβ and the [O III] λλ 5007, 4959 doublet emission lines. The galaxy component (bulge + disc) is shown with a dotted orange line, and is based on Bruzual & Charlot (2003) models. The AGN component, which consists on a full SED template from Temple et al. (2021) for type 1 and narrow emission lines for type 2 is shown with a dashed violet line. The cyan line shows the composed SED. Dust absorption is applied in all cases.

3.5.2 Type 2 spectral energy distributions

For objects flagged as type 2 AGN, typical narrow emission lines from AGN are added to the SED of the bulge component of the host galaxy. The emission lines are modelled with a Gaussian profile, with a total luminosity that was computed using photoionisation models made by Feltre et al. (2016). These models were generated using a standard photoionisation code (CLOUDY, Ferland et al. 2017), and in them, the luminosity of the lines depends on a series of input parameters that describe the physical properties of the region where the narrow lines of AGN are emitted (the narrow-line region, or NLR), and which we explain in detail below:

Gas metallicity: log(O/H) + 12. The gas-phase oxygen abundance is chosen randomly in the range 8.7 < log(O/H) + 12 < 9.3, where we adopted the value log(O/H)_⊙ + 12 = 8.71 for the (gas-phase) solar metallicity, as in Feltre et al. (2016) and Gutkin et al. (2016) for a corresponding dust-to metal mass ratio ξ_d = 0.3 (see below). The fact that we are using only models with solar or super-solar metallicities for the NLR of type 2 AGN is motivated by the fact that these are the models that best sample the region covered by local Seyfert galaxies in standard emissionline diagnostic diagrams (Feltre et al. 2016). These values are also consistent with those found in local AGN from observations (e.g. Peluso et al. 2023).

Ionising spectrum. The ionising spectrum used in the CLOUDY models to represent the accretion disc of the AGN has the shape S_v ∝ ν^α for the wavelength range 0.001 ≤ λ/µm ≤ 0.25, where the power-law index α is an adjustable parameter. We used only models with either α = –1.4 or α = –1.7, randomly chosen, which sample the centre of the range modelled by Feltre et al. (2016); that is, –1.2 to –2.

Ionisation parameter: log U. The ionisation parameter is defined as the dimensionless ratio of the number density of H- ionising photons to that of hydrogen. Using a combination of photoionisation models and high-resolution cosmological zoomin simulations of galaxies, Hirschmann et al. (2017) found that, at fixed stellar mass, U is one of the main physical parameters driving the cosmic evolution of optical-line ratios. To reproduce this effect in our catalogue, the ionisation parameter is selected randomly (evolving with redshift) within the following ranges⁷: $$\{\begin{array}{l}-5<\log U<-3\hfill \\ -4<\log U<-2\hfill \\ -3.5<\log U<-1.5\hfill \end{array}\begin{array}{ll}\text{for}\hfill & 0<z<1\hfill \\ \text{for}\hfill & 1<z<2\hfill \\ \text{for}\hfill & z>2,\hfill \end{array}$$(9)

where these ranges have been derived from Hirschmann et al. (2017) (their Fig. 6, central panel). We note that this is the only NLR parameter in our simulation for which we assumed an evolution with redshift, and therefore the redshift evolution of AGN narrow-line ratios is purely linked to that of log U.

Hydrogen number density, n_H : the volume-averaged hydrogen density of the narrow-line region. Chosen randomly n_H = 10³ or n_H = 10⁴ cm^–3, which are typical gas densities estimated from optical line-doublet analyses of NLRs (see e.g. Osterbrock & Ferland 2006; Binette et al. 2024).

Dust-to-metal mass ratio, ξ_d, which accounts for the depletion of metals onto dust grains in the ionised gas. We used only models with ξ_d = 0.3, which implies assuming that 30 per cent by mass of all heavy elements are in the solid phase (Feltre et al. 2016; Gutkin et al. 2016).

In these models, the intensity of the lines is scaled to the accretion luminosity, L_acc, of the AGN; that is, the luminosity due to the accretion onto the central black hole. Assuming that the bolometric luminosity coming from the AGN, L_bol, is the sum of L_acc and L_X, and that L_bol can be retrieved from the X-ray luminosity with a bolometric correction (L_bol = k_bolL_X), L_acc can be deduced from the following equation: $$L_{\text{acc}}=L_{\text{X}}\left(k_{\text{bol}}-1\right),$$(10)

where again we chose k_bol = 25 for consistency with the rest of the work.

The emission lines of type 2 AGN are characterised by velocity dispersions that are higher than those of SF galaxies (with FWHM lower than a few hundred km s^-1 ), but lower than those of their type 1 counterparts (with FWHM ≳ 1000 km s^–1). To model this, we used the results from Menzel et al. (2016), who studied the spectroscopic properties of a sample 2578 X- ray selected AGN in the redshift range z = [0.02, 5.0]. We modelled the FWHM distribution of the Hβ emission line (for FWHM_Hβ < 1000 km s^–1) shown in their Fig. 6 as a lognormal distribution centred at 355 km s^–1, covering the range FWHM ~ [200, 1000] km s^–1. We assigned randomly the velocity of dispersion of type 2 AGN in our catalogue following this distribution.

Finally, it is worth noting that we did not add an AGN component to the continuum of the host galaxies of type 2 AGN. While at UV and optical wavelengths this approach can be a good approximation for galaxies with strongly attenuated AGN emission, the AGN contribution is expected to dominate in the mid and far IR, even for the most attenuated sources, due to the emission from the dusty torus surrounding the accretion disc. The inclusion of such an AGN component is left to future work.

The full list of narrow-region lines added to the spectra of type 2 AGN is reported in Table 2. We show in Fig. 9 an example of type 2 AGN with these lines. In a similar fashion as we did for type 1 AGN, we quantified the dominance of the type 2 AGN component with respect to the total (AGN + host galaxy) using the ratio of the emission line flux of [O III]_λ5007 from the AGN component with respect to the total. This quantity is shown also in Fig. 9.

4 Validation of the catalogue

In this section we compare the physical properties of the AGN in our catalogue with data that were not used for its calibration, in order to test the validity of its predictions.

Fig. 10 Fraction of AGN at z < 0.33 from our catalogue (purple triangles), selected with log λsBHAR < 1.5 to compare them with the fraction reported in Birchall et al. (2022) (yellow squares). The error bars of the MAMBO AGN fraction show the Poissonian uncertainty.

4.1 AGN fraction

The AGN HGMF that we used to calibrate the fraction of AGN over the total galaxy population at z < 2 (Bongiorno et al. 2016) is defined for z > 0.30 and ℳ > 10^9.5 M_⊙. To test the validity of the extrapolations we did at lower redshifts and stellar masses, we compare our results with different recent works. For example, in Fig. 10 we compare the fraction of AGN at z < 0.33 in different redshift bins with the results from Birchall et al. (2022), who studied 917 X-ray counterparts of SDSS galaxies. To make our sample as similar as possible to that of Birchall et al. (2022), we further selected the MAMBO AGN with log λ_sBHAR < 1.5 and ℳ > 10^8.5 M_⊙. We see that for all redshift bins the fractions are roughly consistent within their uncertainities.

Regarding low-mass galaxies, Latimer et al. (2021) studied a sample of 495 dwarf local galaxies (ℳ ≤ 3 × 10⁹ M_⊙, z ≤ 0.15) observed by eRosita, and found an upper limit of 1.8% for the AGN fraction. In the same redshift and mass bin our catalogue yields a fraction of 1.4 ± 0.2%, which is roughly compatible with these results. On the other hand, Mezcua & Domínguez Sánchez (2024) found a much higher AGN fraction (about 20%) when studying a sample of dwarf local galaxies from the MaNGA survey using integral-field spectroscopy. We note that the exact AGN fraction at these low stellar mass and redshift regimes is still a matter of debate, and therefore the predictions from our catalogue should be taken with caution here.

4.2 X-ray luminosity

In order to validate the X-ray properties of our catalogue we used two main observables, the X-ray luminosity function and the F_X cumulative number counts. These comparisons are relevant since they give hints about the purity and completeness of our catalogue in comparison to other X-ray selected catalogues.

Figure 11 shows the hard X-ray luminosity function (XLF) of our mock catalogue compared with several observational luminosity functions from the literature. We note that the XLF by Miyaji et al. (2015) was used in B16 as an extra observational constraint to determine their HGMF and SARDF, and therefore our catalogue should reproduce it by construction. In Fig. 11 the XLF from Buchner et al. (2015) is scaled by a factor 0.65 in order to remove the contribution of Compton-Thick (CTK) sources, which were not considered in the work of B16, and therefore are not represented in our catalogue. This factor was chosen because Buchner et al. (2015) found a constant fraction of CTK objects over the total AGN population of about 35% independent of z and L_X. We observe in general a good agreement between the XLF from our mock and the observed ones at all redshift bins.

As a further check, in Fig. 12 we show the cumulative number counts of objects above a given X-ray flux. For this, we estimated the X-ray flux from L_X using $$F_{\text{X}}=\frac{L_{\text{X}}}{4\pi D_L^{2}K(z)},$$(11)

where D_L is the luminosity distance and K(z) is a K-correction of the form $$K(z)={(1+z)}^{\text{Γ}-2},$$(12)

where Γ is the slope of the X-ray spectrum. For the K-correction we assumed a photon index of Γ = 1.4, which corresponds to the slope of the cosmic X-ray background, and therefore should represent the full population of both obscured and unobscured objects. By comparing the N(> F_X) from our catalogue in Fig.12 to different data from the literature, we see overall a good agreement over a large range of flux until $F_{\text{X}}\underset{\sim }{\gtrsim }2\times {10}^{-17}$ erg s⁻¹ cm⁻²; that is, two orders of magnitudes fainter than the sample used to calibrate our methodology.

4.3 Narrow emission lines

Diagnostic diagrams are a frequently utilised tool for identifying AGN. By comparing the ratios of various emission lines, we can gain insight into whether star formation, AGN or a composite of both processes dominate in the spectra of a given galaxy. In this section we employ various optical nebular line diagnostics in order to validate the AGN on our catalogue.

One of the diagrams used for this purpose is the BPT diagram (Baldwin et al. 1981), which plots the ratio of the optical lines [O III]_λ5007/Hβ against [N II]_λ6584/Hα. In this diagram, the Kauffmann et al. (2003) and Kewley et al. (2001) lines are commonly used to separate between objects dominated by AGN emission (that fall on the top-right of the diagram), those dominated by star-formation (in the bottom-left) and those that have composite spectra (in the central region).

We show in Fig. 13 an example of such a diagram for the objects of our catalogue, selected with z ≤ 0.8, since, as was shown by Kewley et al. (2013b), this would be the maximum red-shift up to which it is safe to use the BPT diagram to separate SF galaxies from AGN. For clarity of visualisation, we show only galaxies and AGN selected with observed magnitude m_H < 24. The AGN in this figure are colour-coded according to the ratio of the emission line flux of [O III] from the AGN component with respect to the total (AGN + host galaxy). We remind the reader that in our mock, the emission line flux in type 2 AGN is directly proportional to L_X, and therefore this ratio is also correlated to L_X. We also show as comparison the observed line ratios from local galaxies from the SDSS catalogue. We see that, in general, the simulated AGN from our mock fall in the regions that correspond to AGN-dominated or composite objects. While there are some that fall in the region of SF objects, the great majority of them are AGN with dominant [O III] emission from the host galaxy (and most of them have low X-ray luminosities, L_X < 10⁴² erg s⁻¹). In fact, different studies have pointed out that the BPT diagram is biased towards more luminous AGN, missing objects with low X-ray (Birchall et al. 2022) or optical luminosity (Schawinski et al. 2010).

An alternative emission line diagnostic to classify AGN is the mass–excitation (MEx) diagram (Juneau et al. 2011), which plots [O III]_λ5007/Hβ against M and was calibrated to separated SF galaxies from AGN at 0.3 < z < 1. Figure 14 shows the MEx diagram for the lighcone presented in this work, where both galaxies and AGN have been selected with 0.3 ≤ z ≤ 1.0 and $m_{H_{\text{E}}}<24$. The limit in magnitude allows for a better comparison with the work from Juneau et al. (2011), as it removes the low stellar mass tail of our catalogue. We see that all of the AGN that fall in the AGN locus of the diagram have line emission dominated by the AGN component (at least for the [O III] line). However, there are some AGN-dominated sources that are classified as MEx-SF galaxies. Similarly to the previous discussion regarding the BPT diagram, the majority of them have L_X < 10⁴² erg s⁻¹, while the MEx diagram was validated using AGN selected above this threshold, and therefore, it is difficult to draw clear conclusions from this sub-sample.

Fig. 11 Hard X-ray luminosity function of the total population of AGN of our catalogue in different redshift bins, shown by the pruple line. The shaded region represents the Poissonian uncertainty and the dotted horizontal line marks the limiting density from our lightcone (corresponding to 1 object Mpc−3 dex−1). For comparison, we show several observed XLFs: Ueda et al. (2014), Buchner et al. (2015), Miyaji et al. (2015), Aird et al. (2015) and Vito et al. (2017).

4.4 AGN colours and the ultraviolet luminosity function

An important validation for the catalogue if we intend to reproduce the observed AGN population is the luminosity function. In this section we study the redshift evolution of the AGN UV luminosity function (UVLF) at rest-frame wavelength λ = 1450 Å, where the majority of UV rest-frame data on AGN is gathered and where type 1 AGN typically present a peak in their SED (the so-called “big blue bump”).

In Fig. 15, we show the UVLF of the type 1 AGN from the lightcone presented in this work, in different redshift bins, up to z = 6. We have used the absolute rest-frame magnitude from the GALEX FUV filter (see Table 1) as a proxy for M₁₄₅₀. We note that the uncertainty shown is only Poissonian, and therefore constitutes a lower boundary since it does not include other systematic effects like selection effects or completeness level, which would increase the uncertainty.

In Fig. 15, we also compare our UVLF with different works from the literature. Specifically, we show at all redshift bins the QSO UVLF from Manti et al. (2017), who parameterised the LF both as a DPL and a Schechter function using a collection of state-of-the-art measurements from z = 0.5 to z = 6.5, and also from Kulkarni et al. (2019), who used a sample of more than 80 000 colour-selected AGN from redshift z = 0 to 7.5 to parameterise the UVLF as DPL evolving with redshift. Both of these works use the absolute monochromatic AB magnitude at a rest-frame wavelength of 1450 Å to construct the UVLF. From z = 3 to 6 we show also the UVLF from Finkelstein & Bagley (2022), who studied jointly the UVLF of galaxies and QSO and parameterised each population individually with a modified DPL, in order to account for the drop in the faint end of the LF. Additionally, in the figure we show with horizontal dotted line marks the limiting density from our lightcone (corresponding to 1 object/Mpc³/mag) for each redshift bin. On the other hand, the vertical dotted line shows the break magnitude, M_*, at 1500 Å of the galaxy UVLF; that is, the magnitude where the galaxy contribution to the ionising background could be relevant, and indeed the galaxy number density higher than the UV/optically selected QSO one. Following Parsa et al. (2016), Ricci et al. (2017b), this magnitude is calculated as $$M_{*}={(1+z)}^{0.206}\left(-17.793+z^{0.762}\right).$$(13)

By comparing the UVLF constructed from our catalogue with the ones derived directly from observations, we observe a general agreement up to z ≲ 5, for magnitudes brighter than the break magnitude M_*. At fainter magnitudes (M_UV > M_*) we observe a drop in our QSO UVLF, as was expected from the definition of M_* (see above). We note also that at these faint magnitudes there is a big discrepancy also among the observed QSO UVLFs. For example, from z = 3 to 5, the faint end of our UVLF is orders of magnitude lower than that of Manti et al. (2017) of Kulkarni et al. (2019), but agrees quite well with the only-QSO LF of Finkelstein & Bagley (2022).

On the other side, our mock produces an overprediction of the bright end of the UVLF. This can be partially due to the fact we assigned the AGN/galaxy fraction starting from X-ray selected catalogues, which tend to be more complete than UV/optical selections of AGN.

Colour–colour diagrams that use UV to mid-infrared colours can be used to select AGN from a galaxy and AGN sample. Additionally, these selections have the advantage that they can be quickly applied to very large data sets without spectroscopic information. In this section we use such diagrams to validate the colour properties of the AGN in our mock. For this, we checked different colour–colour diagrams that are known to separate AGN from galaxies. For UV/optical bands, these diagrams are able to select mainly luminous type 1 AGN, since type 2 tend to be completely dominated by the host-galaxy emission at these wavelengths.

In Fig. 16, we show the (i – H) vs. (u – z) diagram, which was found in Euclid Collaboration (2024a) to be the best colour selection to separate type 1 AGN from galaxies using Euclid and Rubin/LSST filters. These filters are described in Table 1. The filled contours show the distribution of type 1 AGN from the lightcone here presented, separated into those dominated by the AGN component in the rest-frame FUV (f_AGN⁸ > 0.5) or the galaxy component (f_AGN < 0.5). We also show the distribution of spectroscopically confirmed Type 1 AGN from the Chandra-COSMOS catalogue by Marchesi et al. (2016a) as a comparison. In all cases, the distributions have been cut at z < 3. We observe a good agreement between the Chandra-COSMOS AGN and those with /_agn > 0.5, while both distributions are clearly separated from that of AGN with f_AGN < 0.5. Also, the MAMBO AGN dominated by the AGN component are in agreement with the best colour by Euclid Collaboration (2024a).

Fig. 12 Cumulative number counts of X-ray flux in the hard band. The violet shaded area corresponds to the AGN in the MAMBO catalogue, and it shows the uncertainty of our values estimated as the Poissonian error. We show as comparison the number counts reported in different works. References: Luo et al. (2017), Gilli et al. (2007), Georgakakis et al. (2008), Marchesi et al. (2016b).

Fig. 13 BPT diagram for galaxies and AGN from the MAMBO light-cone selected with z ≤ 0.8 and $m_{H_{\text{E}}}<24$. Galaxies are shown in grey, while AGN are colour-coded from blue to red according to the ratio of the emission line flux of [O III] from the AGN component with respect to thetotal(AGN+hostgalaxy). Objects from theSDSSDR8 areshown in green in the background, and are traditionally classified as SF galaxies if the fall at the left of the Kauffmann et al. (2003), as AGN if they fall at the right of the Kewley et al. (2001) line, or as composite if they fall in between these two lines.

Fig. 14 Mass–excitation (MEx) diagram for galaxies and AGN from the MAMBO lightcone selected with 0.3 ≤ z ≤ 1.0 and $m_{H_{\text{E}}}<24$. Galaxies are shown in grey, while AGN are colour-coded from blue to red according to the ratio of the emission line flux of [O III] from the AGN component with respect to the total (AGN + host galaxy). The solid black lines show the empirical division found in Juneau et al. (2011) to separate AGN (above the line) from SF galaxies (below the line) and composite galaxies (between the two lines).

Fig. 15 UV luminosity function of type 1 AGN from our catalogue (purple points) compared with different literature LFs in different red-shift bins. The purple shaded area shows the uncertainty of our values, estimated as the Poissonian error. In the case of MAMBO, the UV magnitude is computed in the GALEX FUV filter, while for the other cases, it refers to M1500. The solid orange and dashed green lines show the parametric LF from Manti et al. (2017), parameterised as a DPL and a Schechter function, respectively. The dash-dotted red and black lines show the DPL parametric LFs from Kulkarni et al. (2019) and Finkelstein & Bagley (2022), respectively. The vertical dotted line shows the break magnitude at which the galaxy contribution should start dominating the UVLF (Parsa et al. 2016), while the horizontal dashed line marks the limiting density from our lightcone (corresponding to 1 object/Mpc3/mag) for each redshift bin.

Fig. 16 Colour-colour diagram to separate galaxies from type 1 AGN as in Euclid Collaboration (2024a). Filled contours show the distribution of type 1 AGN from the lightcone here presented, separated into those dominated by the AGN component (fAGN > 0.5, orange contours) or the galaxy component (fAGN < 0.5, blue contours). The distribution of spectroscopically confirmed type 1 AGN from the Chandra-COSMOS catalogue Marchesi et al. (2016a) is shown with non-filled red contours. In all cases, the samples have been cut at z < 3. The contour levels represent iso-density lines, corresponding to the 50th, 75th, 90th, and 99th percentiles of the distribution. The grey line shows the best selection criteria found in Euclid Collaboration (2024a) to separate type 1 AGN from galaxies and type 2 AGN.

5 Uses for future surveys

This section aims to present some examples of how this catalogue can be used to make predictions that shall aid in the preparation of future large surveys (Euclid, Rubin/LSST, Moons, etc.). We focus on two examples, namely the expected number densities of galaxies and AGN in a given photometric band, and the spectroscopic selection of AGN through diagnostic diagrams. For this purpose, we choose the Euclid mission as a case study.

5.1 Caveats

First of all, we would like to warn the reader of a series of caveats to be taken into account when using this catalogue to perform this kind of scientific analysis. Most of them have already been discussed in the text in different sections, but we gather them here for clarity:

The DM lightcone out of which the galaxy and AGN catalogue is constructed contains DM haloes up to z = 10. However, the relations we used to construct the galaxy and AGN catalogue are calibrated at lower redshifts. For example, the galaxy SMF is based on observations up to z = 7.5, the AGN HGMF up to z = 2.5, and the accretion rate distribution up to z = 4. Therefore, the extrapolations we did at higher redshifts are to be taken with precaution.

The spectra of type 2 AGN are constructed by simply adding narrow AGN lines to the continuum of the host galaxy, without adding any nuclear (non-stellar) continuum contribution. This represents a strong assumption that might not hold true for all sources. Furthermore, we separated the AGN in our catalogue into two only categories: type 1 and 2. In reality, however, we know that such binary classification does not represent the full population of AGN, and AGN are sometimes classified in intermediate classes (1.5, 1.9, etc.).

Likewise, we did not include the IR component from the torus, which usually dominates the continuum emission of type 1 and 2 AGN for wavelengths redder than 2 µm.

Throughout this paper we assumed a constant bolometric correction k_bol = 25 (L_bol = k_bolL_X) to ensure self-consistency. However, the bolometric correction is known to correlate with L_X , and can have a wide range of values covering orders of magnitude (Duras et al. 2020).

Finally, we did not include stars in our mock. Some stars, such as cataclysmic variables, have multi-wavelength features that can lead to their misclassification as QSO. Cataclysmic variables appear as point-like sources with QSO-like colours, and often exhibit X-ray emission. Therefore, our mock does not account for contamination from such stellar objects, which is instead present in observations.

We plan to keep developing our pipeline in order to tackle most of these effects.

5.2 The Euclid Surveys

The European Space Agency Euclid space telescope (Euclid Collaboration 2024g), successfully launched in July 2023, is equipped with two on-board instruments. The VISible instrument (VIS, Euclid Collaboration 2024c) carries a single broadband optical filter, I_E. The Near-Infrared Spectrometer and Photometer (NISP, Euclid Collaboration 2024e), on the other hand, possesses three NIR photometric filters: Y_E, J_E, and H_E (Table 1). The analysis performed in the following subsections is restricted to the H_E filter. Additionally, the NISP instrument is equipped with four different low resolution NIR grisms (with a spectral resolution of R = 380 for a 0.5 arcsecond diameter source): one blue grism (0.92–1.3 µm), and three identical red grisms (1.25–1.85 µm).

During its 6-year mission, Euclid will perform two main surveys: a Wide Survey (EWS), planned to cover an area of 14679 square degrees on the sky, and a Deep Survey (EDS) of 53 deg² (Euclid Collaboration 2024g). The limiting magnitude for point sources detected with a minimum S/N of 5 is m_H = 24 in the EWS (Euclid Collaboration: Scaramella et al. 2022), and at least two magnitudes deeper in the EDS (m_H = 26). During the EWS, the sky will be observed with only one pass of the red grism, which for the Hα line translates into a redshift coverage of 0.9 ≤ z ≤ 1.8 and a flux limit of F_Hα > 2 × 10⁻¹⁶erg s⁻¹ cm⁻². Instead, for the EDS the redshift coverage of Ha spans to 0.4 ≤ z ≤ 1 .8 with a limiting flux F_Hα > 6 × 10⁻¹⁷erg s⁻¹ cm⁻².

5.3 Number densities

One of the most important predictions that can be done with this type of catalogue is the number densities (number of objects per squared degree) of galaxies and AGN in a given magnitude bin. We show in Fig. 17 the number density in the Euclid H-band for both galaxies and AGN, showing with vertical lines the limiting m_H for the Euclid Wide and Deep surveys. We show also for comparison the number density of objects in the COSMOS catalogue. We observe a good agreement between the galaxies from our catalogue and the ones from COSMOS 2015 (Laigle et al. 2016). For the AGN sample, we used again the Chandra COSMOS Legacy Spectral Survey catalogue (Marchesi et al. 2016b). We show in Fig. 17 the m_H distribution of type 1 and 2 AGN (classified into these 2 categories by optical spectroscopy), selected with minimum X-ray flux F_X = 1.9 × 10⁻¹⁵ erg s⁻¹ cm⁻². After applying the same flux limit to type 1 and 2 MAMBO AGN, we observe a good agreement between the respective type 1 and 2 populations in MAMBO and COSMOS.

Additionally, we show in Table 3 the expected surface densities and total (integrated) numbers of galaxies and AGN with different selections, in the EWS and the EDS. First, we performed a photometric selection of sources detected above the limiting magnitude of each survey. Furthermore, the numbers shown in Table 3 consider only AGN selected with intrinsic hard band luminosity L_X > 10⁴² erg s⁻¹ , to remove possible non-AGN X-ray emitters. We compare these numbers with the ones reported in Euclid Collaboration (2024h), who performed a similar analysis to the one presented in this work, with the aim of forecasting the expected surface densities of AGN in the Euclid Surveys. In Tables 4 and 5 of Euclid Collaboration (2024h), the authors report the surface densities of AGN detectable in the H_E band for the Euclid Wide and Deep Surveys. For the EWS, the reported surface densities are 4.5 × 10³ deg⁻², 6.8 × 10²deg⁻² and 1.7 × 10³ deg⁻² for all AGN, type 1 and type 2, respectively.

The corresponding numbers for the EDS are 2.4 × 10³ deg⁻², 8.8 × 10² deg⁻² and 3.5 × 10³ deg⁻². We observe in general a good agreement between these numbers and the ones reported in Table 3, especially for the type 1 AGN, while the numbers differ more for the type 2.

In Table 3 we also show the surface densities of sources with Hα emission line flux above the detection limit of each survey (see Sect. 5.2), and only considering sources within the redshift window at which this line will be observed at each survey. For this part, we only considered narrow-line emitters; that is, galaxies and type 2 AGN, since the values of the emission line fluxes for type 1 AGN are not included in the current version of the catalogue.

Fig. 17 Number density counts in the mH magnitude. Galaxies and AGN from the MAMBO lightcone are represented by filled green and blue histograms, respectively. AGN selected with LX > 1042 erg s−1 are shown with a dotted cyan line. We also show the density counts for galaxies (dotted pink line) and AGN from the COSMOS catalogue (Laigle et al. 2016; Marchesi et al. 2016b), where AGN have been selected with FX > 1.9 × 10−15 erg s−1 cm−2, and separated into type 1 and 2 based on optical spectroscopic criteria (dashed yellow and red lines). The solid green and orange lines show the number density counts of type 1 and 2 AGN from our mock after applying the same cut in FX as for the COSMOS sample. We show with vertical dashed black lines the limiting magnitude of the EWS and EDS. The lower panel show the fraction of MAMBO AGN as a function of mH. The different lines represent different subpopulations of AGN, as in the upper panel.

5.4 Diagnostic diagrams

In Sect. 4, we used two emission line diagnostic diagrams, namely the BPT and the MEx, as validation tools for our light-cone, by studying them at the redshift ranges at which these diagrams have been calibrated. In this section, we studied these same diagrams applying the redshift, magnitude and flux limits corresponding to the Euclid Wide and Deep surveys.

Two aspects should be noted in this regard: first, both these diagrams have been calibrated at low redshift (z ≲ 1), while at higher redshifts the physical properties (metallicity, density, ionising radiation, etc.) of the regions where the lines are emitted are expected to change with respect to local conditions. For example, as studied by different works, the position of AGN in the BPT strongly depends on the gas metallicity, and below Z ~ 0.5 Z_⊙, AGN start to populate the SF region side of the diagram (Groves et al. 2006; Kewley et al. 2013a; Hirschmann et al. 2019). Therefore, these diagrams must be taken more cautiously when used at these redshifts with real data. Secondly, a precise study of the number of AGN that Euclid could select with these diagrams (assuming they work at high z) would involve using Euclid-like spectra with realistic noise and resolution, and official Euclid pipelines for the extraction of the line fluxes, which is beyond the scope of this study. Besides, given the spectral resolution of Euclid, the [N II] and Hα lines are likely to be blended in real Euclid spectra (Euclid Collaboration 2024f).

We show in Fig. 18 the BPT and MEx diagrams for the Wide and Deep Euclid surveys, applying in each case the corresponding selection in emission line flux and magnitude limit, and redshift range where all relevant lines will be observed. The number density of objects for each case is also given in the figure.

6 Summary

We developed an empirical workflow to generate mock catalogues of galaxies and AGN starting from a DM-only simulation. Following Girelli (2021), we populated the DM haloes with galaxies by means of a stellar-to-halo mass relation, developed using a subhalo abundance matching technique on observed SMFs. Galaxies were also separated into quiescent or SF ones, following the relative ratio of the blue and red populations in observed SMFs.

In this paper, we further populated galaxies with AGN following observed host galaxy AGN mass functions at different redshifts and AGN accretion rate distribution functions, which were derived starting from X-ray samples of AGN at z < 4 (Bongiorno et al. 2016; Aird et al. 2018). Following Merloni et al. (2014), we separated AGN into optically unobscured (type 1) or obscured (type 2) ones and assigned a proper SED to each of them. For this, for type 1 AGN we used the parametric SED model by Temple et al. (2021), which accounts for the continuum emission of type 1 AGN, as well as their broad and narrow emission lines. For type 2, instead, we added narrow lines generated using photoionisation models (Feltre et al. 2016) to the host galaxy stellar continuum.

We tested this workflow by applying it to a 3.14 deg² DM Millennium lightcone up to z = 10. The result is a mock catalogue of galaxies and AGN with realistic physical properties and observables (such as broadband rest-frame and observed magnitudes and spectral features), complete at least up to a magnitude of m_H ~ 28 and down to a stellar mass of ℳ ~ 10^7.5 M_⊙. We obtained good agreement between our mock data and state-of-the-art observations such as published L_X luminosity functions (e.g. Buchner et al. 2015; Aird et al. 2015), number counts in different NIR to optical bands (in this paper we showed results only on the H band for simplicity), colour-colour diagrams (using u, z, i, and H bands), and emission line diagnostic diagrams (BPT and MEx).

Finally, we demonstrated how this catalogue can be used to make forecasts for future large surveys, using Euclid as an example. We computed the expected surface densities of type 1 and 2 AGN detectable with a given Euclid broad filter. We show the results for the HE band, forecasting that Euclid will observe about 8.3 × 10² and 5.5 × 10³ deg⁻² type 1 and 2 AGN, respectively, selected with m_H < 24, and 6.1 × 10² (2.6 × 10³) deg⁻² type 2 AGN with narrow-line H_α emission with flux F_Hα > 2 × 10⁻¹⁶ (6 × 10⁻¹⁷) erg s⁻¹ cm⁻² in the EWS (EDS), finding good agreement with other published forecasts. We also give examples of the Euclid view of narrow-line diagnostic diagrams, which are used to separate local AGN from SF galaxies.

The full workflow is designed to be as computationally efficient as possible so that it can be run on a personal computer. In Appendix B, we give more details on the execution time of the main steps of the method.

We plan to update this workflow in the near future in order to tackle the open issues described in Sect. 5.1. For example, we plan to revise the AGN fraction and the assignment of the X-ray luminosity using more up-to-date accretion rate distributions (e.g. Zou et al. 2024). Also, we are working on the inclusion of an obscured AGN continuum for type 2 sources, as well as the emission from a dusty torus, which dominates the IR emission of type 1 and 2 AGN. Additionally, we plan to investigate the clustering properties of the AGN in our mock catalogues and compare them with observational data as an additional validation step. Although similar studies have been conducted for MAMBO galaxies (Girelli 2021), we leave the corresponding analysis for AGN to future work. Finally, we plan to apply this method starting from a bigger DM lightcone.

Fig. 18 Predictions from our lightcone for the BPT (left panels) and MEx diagrams (right panels) as observed by the Euclid Wide and Deep surveys. For each panel we have plotted the galaxies (grey dots) and AGN (colour-coded from blue to red) corresponding to the specific redshift range, magnitude limit (in H-band) and Hα line flux limit of each survey. The surface density of galaxies and AGN in each diagram are given in the figure. See captions of Figs. 13 and 14 for further details.

Data availability

The catalogue used for the analysis performed in this paper, together with the necessary code to produce it, are available through the GitHub repository https://github.com/xalolo/MAMBO/tree/main.

Acknowledgements

We thank the referee for the careful reading and helpful suggestions that helped to improve the manuscript. We acknowledge the helpful conversations and input from Johan Comparat, Andrea Merloni, Matthew Selwood, James Aird, Matthew J. Temple, Roberto Gilli, Marco Mignoli. XLL acknowledges the “IAC International Scholarships Program”. LP and MB acknowledge ASI and Italian Ministry grant “Premiale MITIC 2017”. MB acknowledges the support from INAF Minigrant 2023. AF acknowledges the support from the project “VLT-MOONS” CRAM 1.05.03.07, INAF Large Grant 2022 “The metal circle: a new sharp view of the baryon cycle up to Cosmic Dawn with the latest generation IFU facilities” and INAF Large Grant 2022 “Dual and binary SMBH in the multi-messenger era”. VA acknowledges the support from the INAF Large Grant “AGN & Euclid: a close entanglement” 01.05.23.01.14. IEL acknowledges support from the Cassini Fellowship program at INAF-OAS and the European Union’s Horizon 2020 research and innovation program under Marie Sklodowska-Curie grant agreement No. 860744 “Big Data Applications for Black Hole Evolution Studies” (BiD4BESt).

Appendix A Rejected methodologies

In this appendix we show different tests that we did before arriving at the final workflow presented in this paper and that motivates some of the choices presented above.

Fig. A.1 Probability of a galaxy to be an AGN as in Fig. 4 but using the original Schechter fit from B16 to assign the AGN fraction at z < 1.15. The biggest difference with respect to the distribution shown in Fig. 4 is in the high density of low-mass AGN at low redshift.

In Sect. 3.2 we noted that, when using the HGMF from Bongiorno et al. (2016), we modified the slope of the Schechter fit at the two lowest redshift bins. This decision was motivated because the original Schechter fit from B16, when extrapolated to ℳ < 10^9.5, predicts a large fraction of low-mass AGN. This can be seen in Fig. A.1

In principle, one could avoid doing the steps presented in Sects. 3.2 and 3.3 and instead derive the AGN fraction and the X-ray distribution in one step, starting from an observed accretion rate distribution. For this, we used the p (log λ_sBHAR | ℳ, z) distributions from A18 at all ℳ and z and statistically sample random values from those distributions in order to assign λ_sBHAR to every galaxy in our catalogue, and later LX using Eq. (1). After this, AGN can be selected as objects above a given threshold in either L_X or λ_sBHAR.

We show in Fig. A.2 the X-ray luminosity function resulting from this methodology, both for AGN selected with λ_sBHAR > 0.01 (for consistency with the definition adopted in A18 and in this work), and also for all the X-ray emitters (without any constrain in λ_sBHAR). It is visible from this figure that in both cases, the AGN population coming from this methodology is underestimated when comparing it to the observed ones, especially for z ≲ 1 and z ≳ 3.

Since the p (log λ_sBHAR | ℳ, z) distributions from A18 are defined up to z ≤ 4, we also explored the possibility of using them to infer L_X at all z, starting from the AGN fraction inferred from B16 (as in Sect. 3.2). As it can be seen in Fig. A.3, the results from this methodology tend to overestimate the XLF with respect to observed ones, especially in the range 0.8 ≲ z ≲ 1.8, and therefore we decided to use the approach explained in Sect. 3.3.

Appendix B Computational time

The most computationally expensive step of the workflow described in this paper is running the modified version of the C++ code EGG. In this step, physical properties and observables are assigned to every galaxy and AGN from the catalogue. The execution time scales mainly with two parameters; namely, the number of sources and the number of filters at which one wishes to compute the rest-frame and observed magnitudes. Running EGG on a 2011 iMac with 12 GB of RAM memory and a 2.8 GHz Quad-Core Intel Core i7 processor, with 10⁶ sources and 10 filters both for rest-frame and observed magnitudes took a total time of 36 min 51 s; therefore, around 2.2 ms per source. In all our tests, the running time scales almost linearly with the number of sources.

Fig. A.2 Hard X-ray luminosity function, using the p (log λsBHAR | ℳ, z) distributions from A18 to infer the AGN fraction and LX distribution at all ℳ and z. The dashed brown line corresponds to all the X-ray emitters, while the dashed green line corresponds to the objects selected with λsBHAR > 0.01. For comparison, the solid purple line shows the XLF from our lightcone using the methodology adopted in this paper, as in Fig. 11. See the caption of Fig. 11 for further details and references to the observed XLFs.

Fig. A.3 Hard X-ray luminosity function, using the AGN fraction inferred from B16 (as in Sect. 3.2), but the p (log λsBHAR | ℳ, z) distributions from A18 to infer LX at all z (dashed orange line). For comparison, the solid purple line shows the XLF from our lightcone using the methodology adopted in this paper, as in Fig. 11. See the caption of Fig. 11 for further details and references to the observed XLFs.

References

All Tables

All Figures

	Fig. 1 Schematic view of the full workflow presented in this paper to create the galaxy and AGN mock catalogue. Yellow boxes represent the result of each step, while blue boxes show the necessary inputs. Steps within the dotted green box take place within the modified version of the public code EGG. A detailed description is given in Sects. 2 and 3.
In the text

	Fig. 2 Distribution of stellar mass and redshift of all the galaxies in the MAMBO lightcone used in this work, separated into the SF (blue) and quiescent (red) populations. The dashed black line shows the total (SF+Q) galaxy population.
In the text

Fig. 3 Stellar mass functions of galaxies and AGN in the three redshift bins studied in B16. The green and purple shaded areas correspond to the galaxies and AGN from the MAMBO lightcone, respectively. We also show the best-fit Schechter function to the HGMF of AGN from Bongiorno et al. (2016) with a dotted black line, and the modified Schechter fit used in this work with a dashed black line. The shaded grey area shows the uncertainty in a. Black dots with error bars show the AGN mass function computed using the Vmax method, derived in B16. The lower pannels show the AGN fraction derived as the ratio of the modified Schechter function and the galaxy mass function from MAMBO (dashed blue line), and the same but using the original Schechter function from B16 (dotted black line).

In the text

Fig. 4 Probability of a galaxy hosting an AGN as a function of redshift. Each point corresponds to a galaxy in the MAMBO lightcone, colour-coded in different ℳ bins indicated in the legend. For comparison, the points with error bars connected by dashed lines show the AGN fraction derived in A18, also in different mass bins. Vertical dotted lines mark the centre of the three redshift bins studied in B16; namely, z = 0.55, 1.15, 2.0.

In the text

	Fig. 5 Specific accretion rate distribution function (SARDF) computed at different values of ℳ using Eq. (3), in the three redshift bins from B16.
In the text

Fig. 6 Distribution of intrinsic hard-band X-ray luminosity (LX) of all the AGN in the lightcone, separated into type 1 (blue) and type 2 (orange) AGN, in different redshift bins. We also show with dashed and solid lines the subpopulations of type 1 and 2 in SF and quiescent galaxies. The vertical dotted line marks LX = 1042 erg s–1, the threshold typically applied to separate X-ray emission from AGN or from other origins.

In the text

Fig. 7 AGN from our lightcone scattered in the LX–ℳ plane, separated into type 1 (blue) and type 2 (orange) AGN. The points in the central panel show the distribution of all the sources, while the solid contours and the histograms from the upper and right panel correspond to sources selected above FX = 1.9 × 10–15 erg s–1 cm–2. We show also the distribution of sources from the Chandra COSMOS Legacy Spectral Survey (Marchesi et al. 2016b), with the same cut in FX, and separated into type 1 (cyan) and type 2 (pink) AGN. The contour levels represen iso-density lines, corresponding to the 50th, 75th, 90th, and 99th percentiles of the distribution. The dash-dotted and dotted black lines mark the locus where λSBHAR = 0.01 and 1, respectively, which assuming a mean bolometric correction of kbol = 25 and a constant mass ratio of black hole to host galaxy of MBH ≈ 0.002ℳ correspond approximately to 1 and 100% of the Eddington limit, respectively.

In the text

Fig. 8 Fraction of optically obscured (type 2) AGN as a function of intrinsic X-ray luminosity. The solid green line shows the relation used to separate AGN (Eq. (4), originally from Merloni et al. 2014) into optically obscured/unobscured, while the dashed purple line shows the actual fraction of type 2 AGN in our catalogue (with the Poissonian uncertainty shown by the shaded area). The other two dashed lines give information about the X-ray obscuration (where type 1 X-ray AGN are defined as objects with logNH/cm–2 < 22 and vice-versa): the dashed orange line shows the fraction of objects classified as type 2 both in optical and the X-ray, while the dashed red line represents the fraction of all the AGN that are type 1 in both bands at the same time.

In the text

Fig. 9 Some examples of SEDs of AGN from our lightcone. We show three representative cases: in the upper panel, we show a type 1 AGN with dominant QSO component in the rest-frame UV band, while in the middle panel, we show a type 1 AGN with similar QSO and host-galaxy contribution for the same band. The lower panel shows a type 2 AGN, with a zoom-in subpanel showing the Hβ and the [O III] λλ 5007, 4959 doublet emission lines. The galaxy component (bulge + disc) is shown with a dotted orange line, and is based on Bruzual & Charlot (2003) models. The AGN component, which consists on a full SED template from Temple et al. (2021) for type 1 and narrow emission lines for type 2 is shown with a dashed violet line. The cyan line shows the composed SED. Dust absorption is applied in all cases.

In the text

	Fig. 10 Fraction of AGN at z < 0.33 from our catalogue (purple triangles), selected with log λsBHAR < 1.5 to compare them with the fraction reported in Birchall et al. (2022) (yellow squares). The error bars of the MAMBO AGN fraction show the Poissonian uncertainty.
In the text

Fig. 11 Hard X-ray luminosity function of the total population of AGN of our catalogue in different redshift bins, shown by the pruple line. The shaded region represents the Poissonian uncertainty and the dotted horizontal line marks the limiting density from our lightcone (corresponding to 1 object Mpc−3 dex−1). For comparison, we show several observed XLFs: Ueda et al. (2014), Buchner et al. (2015), Miyaji et al. (2015), Aird et al. (2015) and Vito et al. (2017).

In the text

	Fig. 12 Cumulative number counts of X-ray flux in the hard band. The violet shaded area corresponds to the AGN in the MAMBO catalogue, and it shows the uncertainty of our values estimated as the Poissonian error. We show as comparison the number counts reported in different works. References: Luo et al. (2017), Gilli et al. (2007), Georgakakis et al. (2008), Marchesi et al. (2016b).
In the text

Fig. 13 BPT diagram for galaxies and AGN from the MAMBO light-cone selected with z ≤ 0.8 and $m_{H_{\text{E}}}<24$. Galaxies are shown in grey, while AGN are colour-coded from blue to red according to the ratio of the emission line flux of [O III] from the AGN component with respect to thetotal(AGN+hostgalaxy). Objects from theSDSSDR8 areshown in green in the background, and are traditionally classified as SF galaxies if the fall at the left of the Kauffmann et al. (2003), as AGN if they fall at the right of the Kewley et al. (2001) line, or as composite if they fall in between these two lines.

In the text

Fig. 14 Mass–excitation (MEx) diagram for galaxies and AGN from the MAMBO lightcone selected with 0.3 ≤ z ≤ 1.0 and $m_{H_{\text{E}}}<24$. Galaxies are shown in grey, while AGN are colour-coded from blue to red according to the ratio of the emission line flux of [O III] from the AGN component with respect to the total (AGN + host galaxy). The solid black lines show the empirical division found in Juneau et al. (2011) to separate AGN (above the line) from SF galaxies (below the line) and composite galaxies (between the two lines).

In the text

Fig. 15 UV luminosity function of type 1 AGN from our catalogue (purple points) compared with different literature LFs in different red-shift bins. The purple shaded area shows the uncertainty of our values, estimated as the Poissonian error. In the case of MAMBO, the UV magnitude is computed in the GALEX FUV filter, while for the other cases, it refers to M1500. The solid orange and dashed green lines show the parametric LF from Manti et al. (2017), parameterised as a DPL and a Schechter function, respectively. The dash-dotted red and black lines show the DPL parametric LFs from Kulkarni et al. (2019) and Finkelstein & Bagley (2022), respectively. The vertical dotted line shows the break magnitude at which the galaxy contribution should start dominating the UVLF (Parsa et al. 2016), while the horizontal dashed line marks the limiting density from our lightcone (corresponding to 1 object/Mpc3/mag) for each redshift bin.

In the text

Fig. 16 Colour-colour diagram to separate galaxies from type 1 AGN as in Euclid Collaboration (2024a). Filled contours show the distribution of type 1 AGN from the lightcone here presented, separated into those dominated by the AGN component (fAGN > 0.5, orange contours) or the galaxy component (fAGN < 0.5, blue contours). The distribution of spectroscopically confirmed type 1 AGN from the Chandra-COSMOS catalogue Marchesi et al. (2016a) is shown with non-filled red contours. In all cases, the samples have been cut at z < 3. The contour levels represent iso-density lines, corresponding to the 50th, 75th, 90th, and 99th percentiles of the distribution. The grey line shows the best selection criteria found in Euclid Collaboration (2024a) to separate type 1 AGN from galaxies and type 2 AGN.

In the text

Fig. 17 Number density counts in the mH magnitude. Galaxies and AGN from the MAMBO lightcone are represented by filled green and blue histograms, respectively. AGN selected with LX > 1042 erg s−1 are shown with a dotted cyan line. We also show the density counts for galaxies (dotted pink line) and AGN from the COSMOS catalogue (Laigle et al. 2016; Marchesi et al. 2016b), where AGN have been selected with FX > 1.9 × 10−15 erg s−1 cm−2, and separated into type 1 and 2 based on optical spectroscopic criteria (dashed yellow and red lines). The solid green and orange lines show the number density counts of type 1 and 2 AGN from our mock after applying the same cut in FX as for the COSMOS sample. We show with vertical dashed black lines the limiting magnitude of the EWS and EDS. The lower panel show the fraction of MAMBO AGN as a function of mH. The different lines represent different subpopulations of AGN, as in the upper panel.

In the text

Fig. 18 Predictions from our lightcone for the BPT (left panels) and MEx diagrams (right panels) as observed by the Euclid Wide and Deep surveys. For each panel we have plotted the galaxies (grey dots) and AGN (colour-coded from blue to red) corresponding to the specific redshift range, magnitude limit (in H-band) and Hα line flux limit of each survey. The surface density of galaxies and AGN in each diagram are given in the figure. See captions of Figs. 13 and 14 for further details.

In the text

	Fig. A.1 Probability of a galaxy to be an AGN as in Fig. 4 but using the original Schechter fit from B16 to assign the AGN fraction at z < 1.15. The biggest difference with respect to the distribution shown in Fig. 4 is in the high density of low-mass AGN at low redshift.
In the text

Fig. A.2 Hard X-ray luminosity function, using the p (log λsBHAR | ℳ, z) distributions from A18 to infer the AGN fraction and LX distribution at all ℳ and z. The dashed brown line corresponds to all the X-ray emitters, while the dashed green line corresponds to the objects selected with λsBHAR > 0.01. For comparison, the solid purple line shows the XLF from our lightcone using the methodology adopted in this paper, as in Fig. 11. See the caption of Fig. 11 for further details and references to the observed XLFs.

In the text

Fig. A.3 Hard X-ray luminosity function, using the AGN fraction inferred from B16 (as in Sect. 3.2), but the p (log λsBHAR | ℳ, z) distributions from A18 to infer LX at all z (dashed orange line). For comparison, the solid purple line shows the XLF from our lightcone using the methodology adopted in this paper, as in Fig. 11. See the caption of Fig. 11 for further details and references to the observed XLFs.

In the text