© The Authors 2025

1. Introduction

Massive black holes (MBHs) with masses of 10⁵ − 10⁹ M_⊙ reside at the centre of most galaxies in our local Universe (Gehren et al. 1984; Kormendy & Richstone 1995; Reines et al. 2011, 2013; Baldassare et al. 2020), including in our own Milky Way. These cosmic giants play a crucial role in the evolution and dynamics of galaxies, as was demonstrated by numerous scaling relations between the properties of galaxies and their MBHs (Ferrarese & Merritt 2000; Tremaine et al. 2002; Kormendy & Ho 2013; McConnell & Ma 2013; Schramm & Silverman 2013; Scott et al. 2013; Shankar 2013; Shankar et al. 2020; Sahu et al. 2019; Greene et al. 2020), such as the M_* − M_BH relation with the stellar mass of the host galaxy. These linked properties indicate that MBHs track the hierarchical formation of galaxies, and they might indicate that MBHs and galaxies mutually affect each other. The details of the processes governing the origin, growth, activity, and feedback of these black holes (BHs) remain highly debated, however.

The general scenario is as follows. MBHs must have originated at high redshift from lighter black hole seeds. Several seed models have been proposed that ranging from “light” to “heavy” masses (for a review, see Latif & Ferrara 2016; Volonteri et al. 2021). At the largest scales, the fate of a seed is determined by the evolution of its host dark matter (DM) halo, which proceeds hierarchically in mass. Within the host galaxy, the growth of the seed proceeds primarily through the accretion of gas (Yu & Tremaine 2002) and in a minor part via the accretion of stars (Rees 1988), with periods of intense accretion activity likely alternating with quiescent phases. Intense accretion emits significant energy in a variety of forms: Radiation is emitted across the electromagnetic spectrum (Fabian 2012), while kinetic and thermal energy can affect the gas temperature, density, and turbulence, and it can in turn prevent new gas from reaching the MBH and temporarily suppress the process (AGN feedback). To a lesser extent (Ricarte et al. 2019), BH mergers following the merger of their host galaxies can also contribute to the growth of MBHs (Dubois et al. 2014a; Kulier et al. 2015).

One of our best probes of the population of MBHs is the luminosity function (LF) of accreting MBHs: quasars and active galactic nuclei (AGNs) more generally. This has been the target of surveys in the X-ray, UV, optical, infrared, and radio bands (see Hopkins et al. 2007; Shen et al. 2020, and references therein). These surveys cover distances of up to z ≲ 7 and revealed that the LF evolves strongly with redshift in normalisation and shape. Interestingly, mid-infrared surveys appeared to differ from these results. Using observations from the Spitzer Space Telescope survey (Werner et al. 2004), Lacy et al. (2015) found an increased LF up to z ∼ 3. The higher-redshift mid-infrared Universe is now being probed by the James Webb Space Telescope (JWST) (Gardner et al. 2006). It has detected dozens of candidate AGNs at high redshift (z ∼ 4 − 11), including some at lower bolometric luminosities than were observed before (L_bol ∼ 10⁴⁴ − 10⁴⁶ erg/s) (e.g. Onoue et al. 2023; Kocevski et al. 2023; Maiolino et al. 2024; Übler et al. 2023; Kokorev et al. 2023; Lyu et al. 2024). The first results for the LF and mass function derived from these observations (Harikane et al. 2023; Maiolino et al. 2024; Greene et al. 2024; Matthee et al. 2024; Taylor et al. 2025) also indicate an increased LF compared to previous expectations based on the compilation of data from mid-IR to X-ray, but dominated at high redshift by X-ray and UV observations (Shen et al. 2020). It has been suggested that the masses and abundance of (candidate) MBHs found by the JWST at very high redshifts might pose a challenge for lighter seed scenarios and might require continuous and intense accretion even for heavier seeds (e.g. Harikane et al. 2023; Maiolino et al. 2024; Kokorev et al. 2023; Greene et al. 2024). At the same time, the models that produce the heaviest seeds appear to struggle to produce enough seeds to explain the abundance of (candidate) MBHs (Regan & Volonteri 2024) (see also Habouzit 2024 for a discussion of the compatibility of different simulations of MBH evolution with these results).

These electromagnetic observations are now being complemented with gravitational-wave (GW) observations because pulsar-timing array (PTA) collaborations are likely close to confirming the detection of a stochastic background consistent with merging MBH binaries with masses ≳10⁸ M_⊙ (Antoniadis et al. 2023a; Agazie et al. 2023a; Tarafdar et al. 2022; Reardon et al. 2023; Xu et al. 2023). It is generally thought that the preliminary result for the PTA background, with its relatively high amplitude, would require MBHs to merge and accrete efficiently (Agazie et al. 2023b; Antoniadis et al. 2024a; Barausse et al. 2023; Izquierdo-Villalba et al. 2024) (although see Goncharov et al. 2024). In the next decade, the Laser Interferometer Space Antenna (LISA) will further extend our understanding of merging MBHs by detecting GWs from binaries with lower masses (10⁴ − 10⁸ M_⊙) and up to very high redshifts (z ∼ 20) (Amaro-Seoane et al. 2017; Colpi et al. 2024).

The astrophysical interpretation of observations is achieved by a comparison to theoretical models for the formation and evolution of MBHs. Cosmological N-body and hydrodynamic simulations have become valuable tools for simulating MBHs through the increased computational power and more accurate sub-grid models for MBH accretion, seeding, feedback, and dynamics (see, e.g. Volonteri et al. 2016; Springel et al. 2017; Kannan et al. 2021; Ni et al. 2022; Bhowmick et al. 2024). In addition, semi-analytical models (SAMs; e.g. Cole et al. 2002; Volonteri et al. 2003; Monaco et al. 2007; Somerville et al. 2008; Benson 2012; Barausse 2012; Ricarte & Natarajan 2018; Bonetti et al. 2019; Dayal et al. 2019; Izquierdo-Villalba et al. 2020; Barausse et al. 2020; Trinca et al. 2022) have been widely used to simulate the cosmic history of MBHs because their computational cost is lower. A computationally cheaper alternative to both SAMs and cosmological simulations are empirical models (see, e.g. Soltan 1982; Small & Blandford 1992; Tucci & Volonteri 2017; Conroy & White 2012; Allevato et al. 2021). Instead of making physical assumptions, these models use observations to empirically (and self-consistently) characterise the evolution of MBHs. Empirical models have grown increasingly more complex, and some (Zhang et al. 2023; Boettner et al. 2025) encompass DM haloes, galaxies, and MBHs and count about 50 parameters. The more complex the model, the greater the number observations that must be used to constrain its parameters, from scaling relations to luminosity and mass functions.

Cosmological simulations and SAMs are too computationally expensive to compare the full range of alternative models and parameter space with observations. This results in a set of discrete models. This is a major obstacle for astrophysical inference within a Bayesian framework. In Toubiana et al. (2021), some of the authors of this paper studied mock LISA data and showed that inference with a finite set of discrete population models could lead to severe biases in the underlying astrophysical parameters. Empirical models, while sufficiently fast for a Bayesian parameter estimation (see, e.g. Zhang et al. 2023; Boettner et al. 2025), focus on statistical and empirical properties and not on the underlying physics.

The goal of this work is to provide a fast parametric approach to infer the physics behind the evolution of MBHs from current observations, and prepare for the LISA mission. We introduce our model POMPOCO, which stands for Parametrisation Of the Massive black hole POpulation for Comparison to Observations. It adopts an approach that is intermediate between SAMs and empirical models by proposing an effective description of the formation and evolution of MBHs within their host haloes. As a first application of POMPOCO, we jointly fit the LF at low and high redshifts and the amplitude of the GW background measured by PTA using a Markov chain Monte Carlo (MCMC) algorithm to determine the parameter sets that are compatible with the desired datasets. We validate our results on the M_* − M_BH relation at z = 0 and at higher redshift, for which we do not explicitly fit. The results agree well with the proposed fit of Greene et al. (2020). Finally, we compute predictions for LISA conditioned on these observations.

This paper is organised as follows. Section 2 describes our parametric model. In Sects. 3 and 4, we construct observables from our model and use them to fit the model parameters to observations. We present our results in Sect. 5 and make predictions for the LISA mission in Sect. 6. We summarise and conclude in Sect. 7. In the remainder of the paper, MBH masses are defined in the source frame.

2. Description of the model

We describe the components of our model below. We also summarise all the parameters and the settings of our model in Table 1. In general, our choice of prescriptions is guided by a combination of physical, observational, and simulation-based considerations, as well as a deliberate effort to maintain flexibility and allow the data to inform the inference. In the current version of our model, we do not track the evolution of the spins of MBHs. We hope to do so in future work.

We adopt a flat ΛCDM cosmology with H₀ = 67.77 km s⁻¹ Mpc⁻¹, Ω_m = 0.3071, Ω_b = 0.04820 and δ_c, 0 = 1.686 for the critical overdensity for spherical collapse at z = 0. Given the tight constraints on cosmological parameters, we expect the uncertainty from cosmology to be small compared to the astrophysical uncertainty.

2.1. Dark matter halo merger tree

We generated DM halo merger trees down to z_max = 20 using the implementation of the extended Press-Schechter formalism (Press & Schechter 1974; Bond et al. 1991; Bower 1991; Lacey & Cole 1993) described in Appendix A of Parkinson et al. (2008). We used updated values given in Benson (2017) for the phenomenological modifications to the original extended Press-Schechter formalism, introduced to better reproduce cosmological simulations. We used a time-varying mass resolution, as proposed in Volonteri et al. (2003): M_h, res = 10⁻³M_h, 0(1 + z)^−3.5, where M_h, 0 is the mass of the halo at z = 0.

The suite of trees was built drawing the halo masses at z = 0 from a log-uniform distribution. Each tree entered the calculation of population properties by reweighting the BHs in that tree by their Press-Schechter weight.

2.2. Massive black hole seeds

We populated leaf haloes – the last haloes along the branches of the merger trees¹ – with seed BHs. We only seeded haloes at z ≥ 10 and with mass above a threshold M_h, seed, with probability f_seed.

The nature and mass spectrum of the original seeds of the observed population of MBHs is currently unconstrained. Proposed seeding mechanisms range from light seeds (below 10³ M_⊙), resulting from the first stars or formed primordially, and heavy seeds (between 10³ M_⊙ and 10⁶ M_⊙) resulting from runaway accretion, mergers or supermassive stars. For a recent review of these mechanisms, which are not mutually exclusive, see Volonteri et al. (2021). We allowed our model to probe the seed distribution by drawing the mass of the seed BH from a truncated log-normal distribution, with mean μ_seed and standard deviation σ_seed. The log-normal distribution is meant to capture a single formation scenario, characterised by a range of seed masses around a central value. At the high end, we limited the mass of the BH seed to 10% of the baryonic mass of the halo M_hΩ_b/(Ω_m − Ω_b). At the low end, we took the minimum seed mass to be 100 M_⊙². We denoted the seeded BH as the primary BH of the DM halo.

2.3. Halo and massive black hole mergers

Following the merger of two haloes, we tracked the evolution of the BHs they contain, distinguishing between major and minor halo mergers, based on the halo mass ratio q_h = M_h, 2/M_h, 1 ≤ 1. Based on the results of our simulations, we set the threshold for major mergers to q_h, major = 0.13, which is similar to the value used in Ricarte & Natarajan (2018).

2.3.1. Formation of black hole binaries

Following major mergers, where q_h ≥ q_h, major, we checked whether either halo hosts a BH. If there was only one BH, this became the primary BH of the newly formed halo. If there was one BH in each, they formed a binary with merger time given by the sum of the dynamical friction (DF) timescale and additional delays (see Sect. 2.3.3). Finally, if there were more BHs, we use the prescription for multiple interactions described below.

In the case of minor mergers, BHs in the lighter halo had to sink into the heavier halo before being counted as primary BH or forming a binary. The BHs of the lighter halo were therefore dubbed as outer BHs of the heavier halo. The sinking time is given by the DF timescale of the satellite halo into the primary one. If the lighter halo contained a binary that merges before the sinking time, we allowed it to merge and assign the merger product to the new halo as an outer BH. We retained (the most massive) half of the outer BHs of the lighter halo, in addition to its potential primary/binary, up to a limit of four outer BHs. We found that this is an acceptable limitation, as small BHs in small haloes play a minor role. Similarly, following major mergers, we retained up to four of the outer BHs of the two haloes, selecting the ones with the shorter sinking times. When outer BHs sunk into the centre of the halo, they either became the primary BH of that halo if it contained none, or formed a binary that merged within t_delay (see Sect. 2.3.3). If the halo already contained a binary, or multiple outer BHs sunk within one time step, we used our prescription for multiple interactions described in Sect. 2.3.2.

2.3.2. Multiple interactions

We used the results of Bonetti et al. (2018a) to handle interactions between more than two BHs. Bonetti et al. (2018a) explored the outcome of three-body interactions between MBHs accounting for relativistic corrections to the Newtonian equations of motion. They found that, on average, there is a probability p_multi = 0.22 (see their Table 2) that the interaction of an MBH binary with a third MBH will lead to a rapid merger (within a few hundred million years), usually between the two most massive MBHs. In the remaining cases, the MBH binary is little affected by the third MBH. Based on these results, we decided with probability p_multi = 0.22 if the BHs undergo a multiple interaction or not. If it occurred, we retained the two most massive BHs and drew the merger time (in yr) from a log-normal distribution with mean 8.4 dex and standard deviation 0.4 dex (see Fig. 7 of Bonetti et al. 2018a). If they did not undergo a multiple interaction, we distinguished between two cases, depending on how the BHs met. If there was a binary and a third BH was brought either by a halo merger or by an outer BH that sunk, we kept the original binary, and its time to merger was unaffected. If there was no binary originally, for instance, if there was a primary BH and two outer BHs sunk within one time step, we just kept the two most massive among them, and formed a binary merging within t_delay (see Sect. 2.3.3). In any case, we did not keep track of the other BHs.

2.3.3. Dynamical friction and other delays

Satellite haloes (and their BHs) sink to the centre of the primary halo via DF, whose typical timescale is set by (Lacey & Cole 1993; Binney & Tremaine 1987)

$$\begin{array}{cc}\hfill t_{\mathrm{DF}}& =0.495\frac{1+q_{\mathrm{h}}}{q_{\mathrm{h}}}\frac{1}{H(z)\sqrt{{\mathrm{\Delta }}_{\mathrm{vir}}}ln(1+q_{\mathrm{h}})},\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{vir}}& =178{\mathrm{\Omega }}_{\mathrm{m}}{(z)}^{0.45},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill q_{\mathrm{h}}& =\frac{M_{\mathrm{h},2}}{M_{\mathrm{h},1}},\hfill \end{array}$$(3)

where we made the same choice for the numerical prefactor as in Volonteri et al. (2003)³, and where H(z) is the Hubble factor at the desired redshift. We neglected any possible variations in the DF timescale due to the initial conditions of the infall, and effects such as tidal stripping and tidal heating of the satellite halo.

Once the satellite halo has sunk into the primary one, the two galaxies they host will merge thanks to additional DF. The BHs at their centre will then form a bound binary, with separation of order of a parsec. The binary will undergo several other processes that affect its orbital decay, including stellar hardening and interactions with nuclear gas (see Colpi 2014, for a review). We modelled these additional delays with a single parameter, t_delay, taken to be the same for all binaries. Our model could be extended by allowing the delay time to be drawn from a parametrised distribution that depends on the physical parameters of the system.

2.3.4. Merger

Once the total time delay had elapsed, binaries merged and formed remnants with mass given by the sum of the masses at the time of the merger. We neglected the mass loss due to the emission of GWs, but we did account for the accretion occurring from binary formation until merger, using the prescriptions described in the section below. We labelled the merger remnant as the primary BH of its host halo.

We did not model kicks on the remnant BH due to the emission of GWs, which could displace the BH from the centre of its host, because they depend sensitively on the BH spins and the details of the orbit at the time of merger (eccentricity, spin inclinations), which we did not track in the current version of POMPOCO. Kicks are expected to have a small impact in high-mass galaxies, because of their large escape velocity. Those are the galaxies that contribute the most to the GW background in the nHz band and the high end of the LF, which are the observables of interest for this paper. Neglecting kicks is therefore expected to have little impact on our results. It might impact the LISA rates, however, because many LISA mergers occur in low-mass galaxies, and could modify the low mass end of the M_* − M_BH relation. We left the inclusion of kicks for future work.

2.4. Accretion

We parametrised the growth of the mass of a BH, m_BH, through accretion by the Eddington ratio, f_Edd,

$$\begin{array}{cc}\hfill {\dot{m}}_{\mathrm{acc}}& =f_{\mathrm{Edd}}{\dot{m}}_{\mathrm{Edd}},\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill {\dot{m}}_{\mathrm{BH}}& =(1-ϵ){\dot{m}}_{\mathrm{acc}},\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill {\dot{m}}_{\mathrm{Edd}}& =\frac{m_{\mathrm{BH}}}{ϵt_{\mathrm{Edd}}}=\frac{L_{\mathrm{Edd}}}{ϵc^2},\hfill \end{array}$$(6)

where ṁ_Edd and L_Edd are the Eddington accretion rate and luminosity, respectively, ϵ is the radiative efficiency, and t_Edd = 450 Myr defines the accretion timescale. In this way, the mass between two time steps increases as

$$\begin{array}{c}\hfill m_{\mathrm{BH}}(t+\mathrm{\Delta }t)=m_{\mathrm{BH}}(t)exp\left(f_{\mathrm{Edd}}\frac{1-ϵ}{ϵ}\frac{\mathrm{\Delta }t}{t_{\mathrm{Edd}}}\right).\end{array}$$(7)

We took the radiative efficiency to be 0.1. Inspired by Ricarte & Natarajan (2018), we considered two accretion modes: burst and steady accretion.

The burst mode was triggered following a major halo merger, q_h ≥ q_h, major. Such mergers are expected to feed the reservoir surrounding MBHs, leading to an episode of intense accretion. When this happened, the most massive BH in the resulting galaxy started accreting for a time t_burst with an Eddington ratio f_Edd, burst. In the burst mode, we did allow for super-Eddington accretion, and drew f_Edd, burst from a power-law distribution with slope (γ_burst − 1) < 0 between 10⁻² and 10. In this work, we considered a common burst time t_burst for all BHs. If the resulting halo underwent another major merger before a time t_burst had elapsed, we restarted the count. The choice of a power-law distribution with negative slope is motivated by the general expectation that super-Eddington accretion should be rare. Aird et al. (2018) found that X-ray observations of AGN activity across redshift and galaxy type suggest a steep power-law at high accretion rates, allowing for rare and short-lived periods of super-Eddington accretion.

Whenever they were not accreting in the burst mode, BHs accreted in the steady mode with an Eddington ratio f_{Edd, steady} drawn from a log-normal distribution with mean μ_steady and variance σ_steady truncated at the high end at 1. When drawing f_{Edd, steady}, we drew jointly t_steady from a log-uniform distribution between 1 and 100 Myr. After a time t_steady had elapsed, we drew again both quantities. This procedure is meant to capture the variability of AGNs. The choice of a log-normal distribution is motivated by Volonteri et al. (2016), Fig. 14, which shows the distribution of accretion rates in AGNs from the Horizon-AGN simulation (Dubois et al. 2014b). Observations also suggest that accretion rates could follow a combination of log-normal and power law distributions, see Aird et al. (2018).

The most massive BHs are observed to be quiescent at low redshift due to gas having been consumed by star formation, as well as affected by supernova and AGN feedback (see Merloni 2004; Hirschmann et al. 2012, 2014, and references therein). To reproduce this, for z ≤ z_cut we shut off accretion for BHs with m_MBH ≥ m_cut(z). We parametrised m_cut(z) as

$$\begin{array}{c}\hfill {log}_{10}{m}_{\mathrm{cut}}(z)={log}_{10}{m}_{\mathrm{cut},0}{(1+z)}^{{\alpha }_{\mathrm{cut}}},\end{array}$$(8)

where log₁₀m_cut, 0 and α_cut are parameters of our model, together with z_cut. We imposed α_cut ≥ 0 so that the mass cut is larger at higher redshift. This approach gives our model the flexibility to capture the anti-hierarchical growth observed in data and simulations (Merloni 2004; Hirschmann et al. 2012, 2014), without enforcing it a priori⁴.

We applied the prescriptions above to the primary and binary BHs of all haloes. Outer BHs do not accrete in our model, as we do not expect them to have a sufficient reservoir to accrete significantly. Both BHs in a binary had the same f_Edd (either burst or steady).

3. Observables

In order to compare our model with observations we generated full MBH population for different choices of parameters. To do so, we applied the following procedure:

• We built N_h merger trees, drawing the mass of the haloes at z = 0 from a log-normal distribution between 10¹¹ and 10¹⁵ M_⊙;

• We populated each merger tree according to the model parameters;

• We saved the properties of the BHs in the merger tree (mass, Eddington ratio, mass of the host halo, etc.) at a set of redshifts.

In this work, we generated N_h = 2000 haloes.

The generation of the DM halo merger tree is the most time-consuming step of our model. Fortunately, none of the parameters affect this step. We could therefore generate a single merger tree and run simulations for many sets of parameters on it. This resulted in a relatively fast model, and allowed us to thoroughly explore the parameter space. As an example, it takes ∼1 hour to generate a merger tree and run the simulation for 500 sets of parameters on it. The 2000 merger trees were run in parallel and the results were combined in post-processing.

We will now describe how we transformed the output of our simulations into quantities that can be compared with observations: the LF, the stellar mass-black hole mass relation, the stochastic GW background, and the merger rate.

3.1. Luminosity function

The luminosity of a BH is given by

$$\begin{array}{c}\hfill L_{\mathrm{BH}}=f_{\mathrm{Edd}}\frac{m_{\mathrm{BH}}}{t_{\mathrm{Edd}}}{c}^2,\end{array}$$(9)

where f_Edd is either the burst mode or steady mode Eddington ratio, depending on the current accretion mode of the BH. We assigned binaries a luminosity equal to the sum of the luminosities of the individual BHs (recall that in our model, BHs in a binary share the same f_Edd). Outer BHs of a halo were not included in the computation of the LF since we do not expect them to accrete significantly.

We denote the LF at redshift z₀ by Φ(L_BH, z₀). It is given by the weighted sum of the contribution from each halo merger tree,

$$\begin{array}{c}\hfill \mathrm{\Phi }(L_{\mathrm{BH}},z_0)={\displaystyle \int W_{\mathrm{PS}}}(M_{\mathrm{h},0})\frac{\mathrm{d}{N}_{M_{\mathrm{h},0}}}{\mathrm{d}{log}_{10}{L}_{\mathrm{BH}}}{|}_{z_0}\mathrm{d}{M}_{\mathrm{h},0},\end{array}$$(10)

where W_PS(M_h, 0) is the Press-Schechter weight for a halo of mass M_h, 0 at z = 0 and $\frac{\mathrm{d}{N}_{M_{\mathrm{h},0}}}{\mathrm{d}{log}_{10}{L}_{\mathrm{BH}}}{|}_{z_0}$ is the number of BHs with luminosity L_BH at redshift z₀ in the merger history of a halo of mass M_h, 0 at z = 0. We computed the integral in a discrete way, by the PS-weighted sum of the individuals histograms of each tree. Binning in log-luminosity, we have for the LF in bin k

$$\begin{array}{c}\hfill \mathrm{\Phi }(L_{\mathrm{BH},k},z_0)=\frac{{\mathcal{D}}_{\mathrm{h}}}{N_{\mathrm{h}}{\mathrm{\Delta }}_{L_{\mathrm{BH}}}}{\displaystyle \sum _i}{W}_{\mathrm{PS}}(M_{\mathrm{h},0,i})\frac{\mathrm{d}{N}_i}{\mathrm{d}{log}_{10}{L}_{\mathrm{BH},k}}{|}_{z_0},\end{array}$$(11)

with Δ_LBH the size of the bin in log₁₀L_BH and 𝒟_h = log₁₀(M_{h, 0, max})−log₁₀(M_{h, 0, min}) coming from the normalisation of the log-uniform distribution used to draw the halo masses. The Poisson error in each bin is obtained by assuming that the histograms resulting from the different simulations are independent Poisson realisations. The variances thus sum up and the error is given by

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Phi }(L_{\mathrm{BH},k},z_0)=\frac{{\mathcal{D}}_{\mathrm{h}}}{N_{\mathrm{h}}{\mathrm{\Delta }}_{L_{\mathrm{BH}}}}\sqrt{{\displaystyle \sum _i}{W}_{\mathrm{PS}}{(M_{\mathrm{h},0,i})}^2\frac{\mathrm{d}{N}_i}{\mathrm{d}{log}_{10}{L}_{\mathrm{BH},k}}{|}_{z_0}}.\end{array}$$(12)

3.2. Relation of the stellar mass to the black hole mass

In our code, we did not directly model the evolution of the stellar mass M_* of a given halo. Instead, we used the halo mass-stellar mass relation of Moster et al. (2010) to convert M_h into M_*, including the scatter in the relation, to account for the spread in the relation. We also accounted for the redshfit evolution of the relation when computing the M_* − M_BH relation at high z.

For a given stellar mass at redshift z₀, we can compute the average black hole mass across haloes at that redshift as

$$\begin{array}{c}\hfill M_{\mathrm{BH}}(M_{\ast },z_0)=\frac{{\displaystyle \int W_{\mathrm{PS}}}(M_{\mathrm{h},0})\frac{\mathrm{d}{N}_{M_{\mathrm{h},0}}}{\mathrm{d}{m}_{\mathrm{BH}}\mathrm{d}{M}_{\ast }}{|}_{z_0}{m}_{\mathrm{BH}}\mathrm{d}{m}_{\mathrm{BH}}\mathrm{d}{M}_{\mathrm{h},0}}{{\displaystyle \int W_{\mathrm{PS}}}(M_{\mathrm{h},0})\frac{\mathrm{d}{N}_{M_{\mathrm{h},\mathrm{z}}}}{\mathrm{d}{M}_{\ast }}{|}_{z_0}\mathrm{d}{M}_{\mathrm{h},0}},\end{array}$$(13)

where $\frac{\mathrm{d}{N}_{M_{\mathrm{h},0}}}{\mathrm{d}{m}_{\mathrm{BH}}\mathrm{d}{M}_{\ast }}{|}_{z_0}$ is the number of haloes at redshift z₀ that contain a BH with mass m_BH and whose stellar mass is M_*, in the merger history of a halo of mass M_h, 0 at z = 0.

In practice, we predefined a fixed binning scheme for the stellar masses, and then associated each halo hosting a BH to the corresponding M_*, k bin based on the M_h − M_* relation of Moster et al. (2010). We have

$$\begin{array}{c}\hfill M_{\mathrm{BH}}(M_{\ast ,k},z_0)=\frac{{\displaystyle {\sum }_i}{W}_{\mathrm{PS}}(M_{\mathrm{h},0,i}){\displaystyle {\sum }_j}{m}_{\mathrm{BH},i,j,k}}{{\displaystyle {\sum }_i}{W}_{\mathrm{PS}}(M_{\mathrm{h},0,i})N_{i,k}},\end{array}$$(14)

where m_{BH, i, j, k} are the masses of the BHs in haloes at z₀ in the i-th simulation whose stellar mass is in the M_*, k bin, and N_i, k is the number of such haloes. In the case of binaries, we used the total mass. Outer BHs are not included in this computation.

We computed the standard deviation of the M_BH − M_* relation as

$$\begin{array}{c}\hfill \mathrm{\Delta }{M}_{\mathrm{BH}}(M_{\ast ,k},z_0)=\sqrt{\frac{{\displaystyle \sum _i}{W}_{\mathrm{PS}}(M_{\mathrm{h},0,i}){\displaystyle \sum _j}{m}_{\mathrm{BH},i,j,k}^2}{{\displaystyle \sum _i}{W}_{\mathrm{PS}}(M_{\mathrm{h},0,i})N_{i,k}}-M_{\mathrm{BH}}{(M_{\ast ,k},z_0)}^2}.\end{array}$$(15)

Technically, there should also be a Poisson error (on N_i), but we neglected it, as we take the bins in stellar mass to be big enough so that the dominant source of error is the variance of BH masses within that bin and not the counting.

3.3. Stochastic gravitational-wave background

We computed the stochastic GW signal in the PTA band of our simulations with the formalism of Phinney (2001) and Sesana et al. (2008). In this work, we assumed binaries to be on quasi-circular orbits. We introduced the chirp mass of a binary with masses m₁ and m₂ as ℳ_c = (m₁³m₂³)/(m₁ + m₂)⁵. Denoting the characteristic GW amplitude at (observer) frequency f as h_c(f), we have

$$\begin{array}{c}\hfill h_c^2(f)=\frac{4G}{\pi f^2{c}^2}{\displaystyle \int \frac{W_{\mathrm{PS}}(M_{\mathrm{h},0})}{1+z}}\frac{\mathrm{d}{N}_{M_{\mathrm{h},0}}}{\mathrm{d}{\mathcal{M}}_c\mathrm{d}z}\frac{\mathrm{d}{E}_{\mathrm{GW}}({\mathcal{M}}_c)}{\mathrm{d}ln{f}_r}\mathrm{d}z\mathrm{d}{\mathcal{M}}_c\mathrm{d}{M}_{\mathrm{h},0},\end{array}$$(16)

where $\frac{\mathrm{d}{N}_{M_{\mathrm{h},0}}}{\mathrm{d}{\mathcal{M}}_c\mathrm{d}z}$ is the number of haloes hosting a binary merger with chirp mass ℳ_c at redshift z in the merger history of a halo of mass M_h, 0 at z = 0, and $\frac{\mathrm{d}{E}_{\mathrm{GW}}({\mathcal{M}}_c)}{\mathrm{d}ln{f}_r}$ is the energy emitted in the source rest-frame at frequency f_r = f(1 + z). For quasi-circular binaries, the latter is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}{E}_{\mathrm{GW}}({\mathcal{M}}_c)}{\mathrm{d}ln{f}_r}=\frac{{(\pi {\mathcal{M}}_{c}G{f}_r)}^{2/3}}{3}{\mathcal{M}}_c.\end{array}$$(17)

The characteristic strain thus has a power-law frequency dependence

$$\begin{array}{c}\hfill h_c^2(f)=\frac{4G^{5/3}}{3{\pi }^{1/3}{c}^2}{f}^{-4/3}{\displaystyle \int W_{\mathrm{PS}}}(M_{\mathrm{h},0})\frac{\mathrm{d}{N}_{M_{\mathrm{h},0}}}{\mathrm{d}{\mathcal{M}}_c\mathrm{d}z}\frac{{\mathcal{M}}_c^{5/3}}{{(1+z)}^{1/3}}\mathrm{d}z\mathrm{d}{\mathcal{M}}_c\mathrm{d}{M}_{\mathrm{h},0}.\end{array}$$(18)

The amplitude is determined here by the distribution of mergers in a given simulation. In practice, we performed a Monte Carlo integration,

$$\begin{array}{c}\hfill h_c^2(f)=\frac{4G^{5/3}{\mathcal{D}}_{\mathrm{h}}}{3{\pi }^{1/3}{c}^2}{f}^{-4/3}{\displaystyle \sum _i}{W}_{\mathrm{PS}}(M_{\mathrm{h},0,i}){\displaystyle \sum _j}\frac{{\mathcal{M}}_{c,i,j}^{5/3}}{{(1+z_{i,j})}^{1/3}},\end{array}$$(19)

where the second sum is performed over all the mergers in the i-th merger tree. The Poisson error is then given by

$$\begin{array}{c}\hfill \mathrm{\Delta }{h}_c^2(f)=\frac{4G^{5/3}{\mathcal{D}}_{\mathrm{h}}}{3{\pi }^{1/3}{c}^2}{f}^{-4/3}\sqrt{{\displaystyle \sum _i}{W}_{\mathrm{PS}}{(M_{\mathrm{h},0,i})}^2{\displaystyle \sum _j}{\left(\frac{{\mathcal{M}}_{c,i,j}^{5/3}}{{(1+z_{i,j})}^{1/3}}\right)}^2}.\end{array}$$(20)

Finally, following Antoniadis et al. (2024a), we define the root-mean-square residual as

$$\begin{array}{c}\hfill {\rho }_c^2(f)=\frac{h_c^2(f)}{12{\pi }^2{f}^3}{\mathrm{\Delta }}_f,\end{array}$$(21)

where Δ_f is the frequency resolution (the inverse of the observation time).

3.4. Merger rate

We introduce m_t = m₁ + m₂ as the total mass of an MBH binary and q_b = m₁/m₂ ≥ 1 its mass ratio. The rate of MBH mergers per unit m_t, q_b and z is computed as

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}t\mathrm{d}{m}_t\mathrm{d}{q}_b\mathrm{d}z}={\displaystyle \int 4\pi }{d}_c{(z)}^{2}c{W}_{\mathrm{PS}}(M_{\mathrm{h},0})\frac{\mathrm{d}{N}_{M_{\mathrm{h},0}}}{\mathrm{d}{m}_t\mathrm{d}{q}_b\mathrm{d}z}\mathrm{d}{M}_{\mathrm{h},0},\end{array}$$(22)

where d_c(z) is the comoving distance at redshift z and $\frac{\mathrm{d}{N}_{M_{\mathrm{h},0}}}{\mathrm{d}{m}_t\mathrm{d}{q}_b\mathrm{d}z}$ is the number of MBH mergers with parameters (m_t, q_b, z) in the merger history of a halo of mass M_h, 0 at z = 0, t is the detector-frame time. Monte Carlo integration gives

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}t\mathrm{d}{m}_t\mathrm{d}{q}_b\mathrm{d}z}=\frac{{\mathcal{D}}_{\mathrm{h}}}{N_{\mathrm{h}}}{\displaystyle \sum _i}4\pi d_c{(z)}^{2}c{W}_{\mathrm{PS}}(M_{\mathrm{h},0,i})\frac{\mathrm{d}{N}_i}{\mathrm{d}{m}_t\mathrm{d}{q}_b\mathrm{d}z}·\end{array}$$(23)

In practice, to generate samples of observed MBH merger events for an observation time T_obs we used the following procedure. We read the (m_t, q_b, z) values of all merger events in our trees, attributed them a Poisson rate

$$\begin{array}{c}\hfill \lambda (m_t,q_b,z)=\frac{4\pi d_c{(z)}^{2}c{W}_{\mathrm{PS}}(M_{\mathrm{h},0,i}){\mathcal{D}}_{\mathrm{h}}}{N_{\mathrm{h}}},\end{array}$$(24)

and drew the number of such events from a Poisson distribution with expected number λ(m_t, q_b, z)T_obs.

The detector-frame rate as a function of a single parameter, for example m_t, was obtained by marginalising over the remaining parameters, and the total rate was obtained by marginalising over all parameters (m_t, q_b, z)

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}t}={\displaystyle \int 4\pi }{d}_c{(z)}^{2}c{W}_{\mathrm{PS}}(M_{\mathrm{h},0})\frac{\mathrm{d}{N}_{M_{\mathrm{h},0}}}{\mathrm{d}{m}_t\mathrm{d}{q}_b\mathrm{d}z}\mathrm{d}{M}_{\mathrm{h},0}\mathrm{d}{m}_t\mathrm{d}{q}_b\mathrm{d}z,\end{array}$$(25)

see also Eq. (2) of Hartwig (2019). Finally, Monte Carlo integration gives

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}t}={\displaystyle \sum _{(m_t,q_b,z)}}\lambda (m_t,q_b,z).\end{array}$$(26)

4. Fitting to observations

As a first application of POMPOCO, we fitted jointly the LF across redshifts and the GW stochastic background.

4.1. Observations

As the observed LF, we considered a combination of the fits from Lacy et al. (2015) (fit: ‘All’) and Shen et al. (2020) (fit: ‘global A’), as described in the next section (Fig. 1). The fit by Shen et al. (2020) is based on a wide range of observations in the rest-frame mid-IR, B band, UV, soft and hard X-ray going up to z ∼ 7. The constraint by Lacy et al. (2015) is only based on the mid-IR Spitzer Space Telescope survey with observations up to z ∼ 3. It is also included in the compilation of Shen et al. (2020), although the latter is mainly driven by the more abundant B, UV, and X-ray data. The mid-IR LF fit in Lacy et al. (2015), if taken alone, suggests more AGN at a given bolometric luminosity. We further extrapolated the Lacy et al. (2015) fit to higher redshifts and lower luminosities than covered by their data.

Fig. 1. Evolution of the luminosity function. In green, we show the prediction of our model when fitting for the LF itself at redshifts 0.25, 0.5, 3 and 6, as well as for the GW background measured by EPTA, displaying the median (green line) and 90% confidence interval (green band). This should be compared with the fits to the observed LF from Lacy et al. (2015) (pink) and Shen et al. (2020) (purple) and the region in between (light pink), where we allow the observed LF to lie in our likelihood. At redshifts z = 4, 5 we also plot constraints from JWST (squares) from Matthee et al. (2024). The lowest luminosity point is shown in grey and not in black, because it is likely affected by incompleteness (see the discussion in Matthee et al. 2024). At redshift z = 7, we plot a lower bound obtained by Greene et al. (2024) using JWST results.

Although we did not use them to constrain our model, we show JWST results at z = 4, 5 taken from Matthee et al. (2024) (as converted to bolometric luminosity by Habouzit 2024) in Fig. 1. Data from Matthee et al. (2024) was broad-line selected from the FRESCO (Oesch et al. 2023) and EIGER surveys, and span z = 4.2 − 5.5. FRESCO and EIGER are spectroscopic slitless surveys with uniform selection in their fields, providing a complete sample with clear volume estimates. These results also hint towards a higher LF at high redshift than expected from the fit of Shen et al. (2020), in particular when taking into account that JWST is likely detecting only a subpopulation of all AGNs, and the obscured fraction could be as high as 80% (Gilli et al. 2022). We differentiate the lowest luminosity data point by plotting it in grey, as it is likely affected by incompleteness, see discussion in Matthee et al. (2024).

Other JWST observations (Harikane et al. 2023; Maiolino et al. 2024; Scholtz et al. 2025; Taylor et al. 2025) suggest an even larger population of faint AGN, although sample selections are not as clean as in slitless surveys. In Fig. 1, we show a lower bound to the LF at z = 7 obtained by Greene et al. (2024) using continuum-based selections. These results point in the direction of an even higher LF, but it should be noted that this type of selection may include also non-AGN sources (Pérez-González et al. 2024; Akins et al. 2024).

For the stochastic GW background, we used the results of the second data release from EPTA (Antoniadis et al. 2023a,b,c, 2024a,b), which are consistent with those of other PTAs (Agazie et al. 2024), see Fig. 2.

Fig. 2. Stochastic GW background predicted by our model (green) and the free spectrum measured by EPTA (blue). The 90% confidence intervals centred at the median for the amplitude of the GW background in the first and last bins, i.e., at frequencies 1/(10 yr) and 1/(1 yr), are $4.8_{-2.7}^{+5.1}\times {10}^{-15}$ and $1.0_{-0.6}^{+1.1}\times {10}^{-15}$, respectively.

In this work, we did not fit for the M_* − M_BH relation, but we perform a cross validation by comparing our results at z = 0 with the fit from Greene et al. (2020) (‘All’, including scatter), which is based on a large compilation of galaxies. At higher redshifts, we compare to the fits from Pacucci et al. (2023), derived from JWST data at redshifts z = 4 to 7.

4.2. Computation of the likelihood and sampling

Given the current uncertainty on the LF, we wished to explore a scenario where the true LF lies somewhere between the fits proposed in Shen et al. (2020) and Lacy et al. (2015). We denote Φ_Lacy(L_BH, z₀) and Φ_Shen(L_BH, z₀) the mean of the fits from Lacy et al. (2015) and Shen et al. (2020) respectively, and σ_Lacy(L_BH, z₀) and σ_Shen(L_BH, z₀) their standard deviations (in log space). For a given redshift, we defined a probability distribution p_LF(log₁₀Φ, L_BH, z₀) that, in each luminosity bin, is flat between min(log₁₀Φ_Lacy − σ_Lacy, log₁₀Φ_Shen − σ_Shen) and max(log₁₀Φ_Lacy + σ_Lacy, log₁₀Φ_Shen + σ_Shen) and beyond these limits falls as a Gaussian with standard deviation σ_Lacy or σ_Shen accordingly.

Our likelihood for the LF of a given set of parameters Λ, ℒ_LF(Λ, L_BH, z₀), was obtained by including the Poisson error on the LF computed from simulations so that

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{LF}}(\mathrm{\Lambda },L_{\mathrm{BH}},z_0)={\displaystyle \int p}(\mathrm{\Phi }|\mathrm{\Lambda })p_{\mathrm{LF}}({log}_{10}\mathrm{\Phi },L_{\mathrm{BH}},z_0)\mathrm{d}\mathrm{\Phi },\end{array}$$(27)

where p(Φ|Λ) is a normal distribution with mean given by Eq. (11) and standard deviation Eq. (12). In practice, we performed a Monte Carlo integration, drawing samples for Φ from this normal distribution, truncating to positive values of Φ⁵. The total log-likelihood at z₀ was obtained by summing over the luminosity bins in a z₀-dependent range:

$$\begin{array}{c}\hfill ln{\mathcal{L}}_{\mathrm{LF}}(\mathrm{\Lambda },z_0)={\displaystyle \sum _{L_{\mathrm{BH},\min }\le L_{\mathrm{BH}}\le L_{\mathrm{BH},\max }}}ln{\mathcal{L}}_{\mathrm{LF}}(\mathrm{\Lambda },L_{\mathrm{BH}},z_0){\mathrm{\Delta }}_{L_{\mathrm{BH}}}.\end{array}$$(28)

Finally, the LF log-likelihood was obtained by summing over the different redshifts we wished to fit:

$$\begin{array}{c}\hfill ln{\mathcal{L}}_{\mathrm{LF}}(\mathrm{\Lambda })={\displaystyle \sum _{z_0}}ln{\mathcal{L}}_{\mathrm{LF}}(\mathrm{\Lambda },z_0).\end{array}$$(29)

Similarly, we used the EPTA posteriors (Antoniadis et al. 2023a,b,c, 2024a,b) on log₁₀ρ_c, denoted by p_PTA(log₁₀ρ_c), to define the PTA likelihood in a given frequency bin f₀ as

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{PTA}}(\mathrm{\Lambda },f_0)={\displaystyle \int p}({log}_{10}(h_c)|\mathrm{\Lambda })p_{\mathrm{PTA}}({log}_{10}{\rho }_c)\mathrm{d}{log}_{10}\rho .\end{array}$$(30)

We took p(log₁₀(h_c)|Λ) to be a normal distribution and the mean and error were computed from Eqs. (19)–(21)⁶. We performed a Monte Carlo integration using the samples of p_PTA(log₁₀ρ_c) taht were provided with the EPTA data release.

In this work, as a fiducial result, we fitted the LF at z = 0.25, 0.5, 3 and 6. For z₀ = 0.25, we fited in the range of luminosities [10¹⁰, 10^11.5] L_⊙, for 0.5, in the range [10¹⁰, 10¹²] L_⊙, while for z₀ = 3 and 6, we fitted in the range [10¹², 10¹⁴] L_⊙ and [10^11.5, 10^13.5] L_⊙ respectively.

Since we assumed MBH binaries to be on circular orbits, the slope of the GW background is constrained to −2/3, see Eq. (18). Fitting our model on the whole GW frequency range would therefore give the EPTA data a too large relative weight in the total fit without giving it the flexibility to better accommodate the GW spectrum as a whole. To facilitate comparison with previous works (Antoniadis et al. 2024a; Barausse et al. 2023) and because the way in which we estimate the PTA spectrum (assuming a normal distribution) is more accurate at low frequencies (Lamb & Taylor 2024), we focused on the frequency bin at f₀ = 1/(10 yr).

Finally, we assigned zero likelihood to simulations that produce MBHs at z = 0 with masses exceeding 10¹² M_⊙, as these are unrealistically massive. In our model, such MBHs can form from heavy seeds, ∼10⁶ M_⊙, which later accrete in burst mode with Eddington ratios reaching up to 10. By z ∼ 7, these MBHs typically have masses around 10⁸ − 10⁹ M_⊙, and since we limit accretion only at lower redshifts, there is no mechanism to prevent them from continuing to grow, potentially reaching masses incompatible with observations. A more physically motivated accretion model – based on the available mass reservoir, rather than on Eddington ratios and accretion times – would likely prevent the formation of these excessively large MBHs. However, since these are extreme outliers in our simulations (i.e., they appear in very few of the simulations that survive fitting to the LF and/or GW spectrum and in small numbers), we adopted this simple approach for now, leaving the implementation of a more refined accretion description for future work.

The priors used on the parameters of POMPOCO are detailed in Table 1. We explored the parameter space with an MCMC using the Eryn sampler (Karnesis et al. 2023). Since it is computationally cheap to run serially several sets of parameters for a given merger tree, we used many walkers (∼500) and ran the sampler only for a few steps (a few hundred).

5. Combined fit to the luminosity function and the gravitational-wave background

Our results for the LF, the GW background and the M_* − M_BH relation at z = 0 are shown in Figs. 1, 2 and 3 respectively. In each case, we present in green our model’s median (solid line) and 90% confidence interval (shaded region), as well as the observational data.

Fig. 3. Relation of the stellar mass and MBH mass at redshift z = 0. We compare our results, in green, with the fit of the local relation from Greene et al. (2020) (‘all’), in blue, and the data also compiled in Greene et al. (2020), in black. Green points show a (random) subsample of MBHs produced in our simulations.

In Fig. 1, the light pink area denotes the region between the fits from Lacy et al. (2015) (pink) and Shen et al. (2020) (purple) where we allow the observational LF to reside. The recovered LF is remarkably consistent with observations at all redshifts – including those for which we do not explicitly fit. Note that our LF vanishes at high luminosities because of the finite number of DM haloes we simulate, which does not allow us to generate the rare MBHs behind the very bright end of the LF.

The GW spectrum (Fig. 2) we recover aligns well with the EPTA data (in blue) though our median lies slightly lower. This suggests our results might be better compatible with the reduced amplitude of Goncharov et al. (2024), obtained by adopting an observationally driven model for the pulsar noise. To facilitate comparison to Agazie et al. (2023b), we report the median of the GW background at two reference frequencies, together with the errors defining a 90% confidence region:

• $h_c[1/(10\mathrm{yr})]=4.8_{-2.7}^{+5.1}\times {10}^{-15}$,

• $h_c[1/(1\mathrm{yr})]=1.0_{-0.6}^{+1.1}\times {10}^{-15}$.

As a comparison, if, instead of performing the fiducial fit to the LF and PTA described in Sect. 4.2, we use POMPOCO to fit only for the LF of Shen et al. (2020) across redshifts, we obtain the following estimates:

• $h_c[1/(10\mathrm{yr})]=3.9_{-2.5}^{+7.0}\times {10}^{-16}$,

• $h_c[1/(1\mathrm{yr})]=8.0_{-5.1}^{+14.3}\times {10}^{-17}$.

The median values in the latter case are more than one order of magnitude below the ones reported above. Within the model itself, there is not so much the possibility of having a much lower LF at high z and producing a high GW background. In physical terms, a higher LF means higher density of BHs and so this is why the model can match the EPTA results. We recall that we assume MBH binaries to be circular, therefore the slope of the GW background spectrum is fixed to be −2/3 (see Eq. (18)). Including eccentricity and/or interactions with the environment would allow for a more flexible and shallower spectrum (Enoki & Nagashima 2007; Bonetti et al. 2018b; Antoniadis et al. 2024a; Ellis et al. 2024a). We plan to explore these possibilities in future work.

The resulting M_* − M_BH relation at z = 0 (Fig. 3) is in good agreement with the fit from Greene et al. (2020) (in blue), particularly in the region where most of the observational data (black points) are concentrated. A random sub-sample of individual MBHs from our simulations (green dots) shows that our model is capable of producing the most massive BHs observed in galaxies with M_* ∼ 10¹¹ M_⊙, while also producing MBHs of 10⁷ − 10^8.5 M_⊙ in galaxies with M_* ≥ 10¹¹ M_⊙. While these MBHs are clear outliers, our model is able to reproduce them, a feature that some other models find challenging, as discussed in Habouzit et al. (2021).

Our model predicts a break in the M_* − M_BH relation for galaxies with mass M_* < 10^10.5 M_⊙. This originates from the break in the M_h − M_* relation. The relation predicted directly by our model, M_h − M_BH, does not show such a break at small MBH masses (see Fig. A.1 in the Appendix). This break is rather introduced by our use of the M_h − M_* relation of Moster et al. (2010), which shows one precisely at M_h ∼ 10¹² M_⊙, corresponding to M_* ∼ 10^10.5 M_⊙. Such a break is also observed in some cosmological simulations (Habouzit et al. 2021), and is imputed to the transition between SN feedback, suppressing MBH accretion in low-mass galaxies and AGN feedback, which leads to self-regulation of MBH accretion in high-mass galaxies and quenches star formation. For M_* < 10¹⁰ M_⊙, the median MBH mass is of order 10⁷ M_⊙ and lies above observations, with a shallower slope. We note, however, that there are only a few observations with M_* < 10¹⁰ M_⊙, so that the fit of Greene et al. (2020) is dominated by massive galaxies, and might not be accurate for low-mass galaxies. Ultimately, most observations in low-mass galaxies are within our 90% confidence region, and even those outside can actually be reproduced as outliers of our distribution (green points). In any case, we have checked that this potential overestimate of MBHs with mass 10⁷ M_⊙ does not impact our results on the PTA background, which is fully dominated by MBHs with mass > 10⁸ M_⊙.

At the massive galaxies end, our M_* − M_BH relation shows a break at M_* ≳ 10¹² M_⊙ and slightly underestimates the mass of MBHs with respect to the observational fit. We note, however, that there is only one data point in this region (M_* ≳ 10¹² M_⊙), so the fit might not be valid for these galaxies. As discussed in Appendix A, we have verified that this potential underestimation of the MBH masses should have a small impact on the stochastic GW background. In short, the reason is that this affects MBHs in more massive haloes, which are heavily suppressed by the Press-Schechter mass function, and therefore contribute little to the background.

We show the M_* − M_BH relation at higher redshift in Fig. 4. For z = 4, 5 and 6, we plot the AGN candidates from JWST observations in the redshift range z = 4 to 7 reported by Harikane et al. (2023), Übler et al. (2023), Maiolino et al. (2024), Li et al. (2025a). Since the M_h − M_* relation from Moster et al. (2010) is calibrated for M_* ≳ 10⁸ M_⊙, we restrict our plot to this stellar mass range. Additionally, the redshift evolution of their fit is only calibrated up to z ∼ 3.5, so our results at higher redshifts are based on extrapolation. In order to mimic selection effects in this dataset, we also show in dark green our relation when restricting to MBHs with L_BH ≥ 10¹¹ L_⊙. Most data points are within the range predicted by our model when imposing this cut. Note that, when computing the M_* − M_BH relation predicted by POMPOCO, we also count MBHs in a binary as a single BH using the total mass of the binary, as described in Sect. 3.2, because it is usually challenging for observations to distinguish a single MBH from a binary at close separations, although some candidates have beed identified (Maiolino et al. 2024). We also show the fit to the local relation of Greene et al. (2020), and the fit of Pacucci et al. (2023), which is based on the same high redshift data points shown in Fig. 4 (except for the latest from Li et al. 2025a). We caution the reader that the underlying intrinsic M_* − M_BH relation at high redshift, once selection biases are accounted for, is still debated (Pacucci et al. 2023; Li et al. 2025b), as we do not know yet if a large population of undetected MBHs with lower mass ratio between M_BH and M_* exists.

Fig. 4. Relation of the MBH mass and stellar mass at redshift z = 3, 4, 5, and 6. The darker green regions represent our results when restricting to MBHs with luminosity larger than 1011 L⊙, in order to mimic selection effects at high redshift. For z = 4, 5 and 6, we show JWST observations in the redshift range z = 4 to 7 (Maiolino et al. 2024; Übler et al. 2023; Harikane et al. 2023; Li et al. 2025a), some of which were used to obtain the fit of Pacucci et al. (2023) (light orange, excluding the more recent data of Li et al. 2025a). For comparison, we also show the fit to the observed local relation (light blue) of Greene et al. (2020).

Overall, our model predicts a population in the same locus as the observations, although the slope is shallower than the fits of Greene et al. (2020) and Pacucci et al. (2023). Indeed, our model predicts a sizeable population of MBHs with mass almost comparable with the host galaxy mass in low-mass galaxies. We leave a deeper investigation of the possible causes for this discrepancy to future work. Conversely, the model does not include many very high-mass BHs (> 10⁹ M_⊙): this is related again to our simulated merger trees, which do not include the massive and rare haloes hosting the MBHs powering quasars at very bright end of the LF. At z = 3 and 4, we can see that our model starts to build up to the steeper local M_* − M_BH relation of Greene et al. (2020).

Overall, our analysis demonstrates that observations of the LF and of the (likely) GW background observed by PTAs can find a consistent description in our physical model, provided we allow for a more abundant population of AGN compared to global fits of the AGN LF based mainly on optical/UV/X-ray, pre-JWST data.

We show the posterior distribution of the parameters of our model in Figs. 5 and 6. For the sake of readability, the first figure shows the parameters controlling the seeding and merger delays, and the second figure those controlling accretion. We do not observe strong correlations between these two groups of parameters.

Fig. 5. Posterior on the model parameters related to seeding and BH mergers when fitting for the LF and GW background. We show the mean and standard deviation of the log-normal distribution of seed BH masses, μseed and σseed; the minimum mass of haloes seeded Mh, seed and the seeding probability fseed; the delay of binary BH mergers (in addition to halo dynamical friction), tdelay.

Fig. 6. Posterior on the model parameters related to accretion when fitting for the LF and the GW background. We show the slope of the power-law distribution γburst and the duration tburst of burst accretion; the mean and standard deviation of the log-normal distribution of steady accretion rates, μsteady and σsteady; the redshift zcut below which accretion is shut off for heavier MBHs; the slope of the increase of the shut-off mass, mcut, 0 at z = 0, with redshift, αcut.

We obtain a relatively good constraint on the properties of haloes that are seeded, coming from the abundance of BHs in the Universe at different redshifts (i.e., the normalisation of the LF). Our model suggests that ∼10% of haloes with M_h ≳ 10⁷ M_⊙ at z ∼ 20 should host an MBH seed⁷. The range of values favoured for M_h, seed is more compatible with the prediction of relatively heavy seed models for the formation of seeds (see Regan & Volonteri 2024, for an overview and discussion). We obtain a broad posterior on the mean seed mass μ_seed, but with a clear preference for μ_seed ≳ 10³ M_⊙, in better agreement with the predictions of non-extreme heavy seed models (Volonteri et al. 2021). The tendency towards large values of σ_seed suggests that we tend to prefer broad initial distributions for the mass seed. This is likely due to the need to accommodate observations spanning a broad range of MBH masses, together with the restrictive shape of the seed mass distribution currently adopted in POMPOCO (a single Gaussian). Allowing for a combination of seeding channels (e.g. each having its own μ_seed and σ_seed) might improve this, but at the cost of increasing the parameter space. We expect that LISA will allow us to better identify the contribution of different seeding channels compared to current observations (Gair et al. 2011; Sesana et al. 2011; Klein et al. 2016; Toubiana et al. 2021).

We find large delays in MBH mergers to be strongly disfavoured, with t_delay ≲ 1 Gyr at 90% confidence. We recall this is an additional delay to the DF timescale following halo mergers and that we have considered a single value for all MBHs, without a mass or redshift dependence. Our findings are in line with those of Antoniadis et al. (2024a), Barausse et al. (2023), Agazie et al. (2023b), which found that PTA observations favour MBH binaries merging efficiently following galaxy mergers.

Moving on to the accretion parameters, the posterior on γ_burst suggests a log-flat distribution for f_Edd in burst mode following major halo mergers, favouring the occurrence of super-Eddington accretion episodes. These burst episodes are estimated to last a few dozen million years, in agreement with the findings from Hopkins et al. (2006), Capelo et al. (2015). On average, burst accretion episodes increase the MBH mass by a factor of 1.2, rising to 2.8 in the case of super-Eddington episodes. The parameters of steady-mode accretion, in turn, are poorly constrained individually. The anti-correlation between μ_steady and σ_steady means, however, that ∼10% of the MBHs that accrete in this mode have f_Edd > 10⁻² and ∼4% have f_Edd > 10⁻¹. When we label these as AGN, this suggests that the AGN fraction is in the range 4 − 10% at ‘cosmic noon’, which in our model corresponds to z ∼ z_cut. This is in good agreement with observations in different parts of the electromagnetic spectrum, which suggest an AGN fraction at z ∼ 1 − 3 in the range 1 − 10% (Aird et al. 2012; Trump et al. 2013; Donley et al. 2012). The strong correlations between γ_burst and t_burst and between μ_steady and σ_steady point to the fact that a more physically driven approach to describe accretion, for instance based on a reservoir tracking the amount of mass available to accrete, might be possible, allowing to reduce the number of parameters in our model. We leave this exploration for future work.

Finally, our results do support the anti-hierarchical hypothesis for the growth of MBHs, as we observe a clear preference for accretion to shut off below z_cut ∼ 2.5 with m_cut, 0 ∼ 6 × 10⁶ M_⊙. Although not very well constrained, the evolution parameter for the cut-off mass with redshift, α_cut, does seem to favour strictly positive values, leading to the estimate that accretion is suppressed at z_cut ∼ 2.5 for MBHs with mass ≳10⁸ M_⊙. The inferred value of z_cut is in good agreement with observations and more complete numerical simulations, which estimate the peak of accretion to be in the range z = 2 − 3 (Merloni 2004; Hirschmann et al. 2012, 2014). Using our M_* − M_BH relation, we see that the preferred value of m_cut, 0 corresponds to a stellar mass of ∼8 × 10¹⁰ M_⊙, in good agreement with where the break in the M_h − M_* relation of Moster et al. (2010) lies, for M_h ∼ 10¹² M_⊙. This break has been interpreted as being due to the combined effect of the feedback from supernovae (in smaller haloes) and AGNs (in larger haloes) (Croton et al. 2006; Shankar et al. 2006; Bower et al. 2006; Somerville et al. 2008).

The fact that the inferred values of the parameters entering POMPOCO are overall in good agreement with findings from observations and/or more complete numerical simulations reinforces our confidence that our model is able to capture the key features driving the formation and evolution of MBHs and that our conclusions on the compatibility of the LF and the GW background are robust.

6. Merger rates and LISA predictions

Next, we turn to the predictions for merger rates and LISA observations based on our fits to the LF and the GW stochastic background. We compute the yearly rate as a function of MBH binary parameters, as well as the total rate (per year), following the procedure described in Sect. 3.4. To estimate the LISA detection rate, we restrict to events with signal-to-noise ratio (S/N) above 10, which we assume to be the threshold for detection in a 4-year mission. The S/N of the sources is computed assuming the “Proposal” noise curve (Amaro-Seoane et al. 2017) and using the IMRPhenomXHM model (García-Quirós et al. 2020) for the waveform. Since we do not track the evolution of spins in our simulations, for simplicity, we attribute a spin 0 to both BHs in a binary when computing S/Ns. We do not expect our results for the expected number of detections to be much affected by this choice, since spins, at the level of the whole population, introduce a second order correction to the detectability of the sources. For each event, we draw the sky location, the phase at coalescence, the polarisation and the inclination angle isotropically.

Figure 7 shows the posterior distribution of the merger rate per year, while Figs. 8 and 9 present this rate as a function of the total binary mass, m_t, and redshift, z, respectively. The thick lines represent the median values, with the shaded regions indicating the 90% confidence intervals, and the thin lines illustrating individual realisations of our model. Green curves show the intrinsic rate, while the light green ones show the rate of events detectable with LISA. Our predictions for LISA cover a wide range, with the expected number of merger events varying from a few to thousands per year.

Fig. 7. Prediction of the yearly rate of MBH mergers after fitting for the LF and the GW background measured by EPTA. The solid curve shows the intrinsic rate, the dotted one the rate of detetable events, assuming an S/N threshold of 10. For comparison, we also report the predictions for the intrinsic rate of other models from the literature (see text for description). For Q3-nod, Q3-d and PopIII-d, the shaded areas show the range reported in reference Barausse et al. (2023) (the lower bound corresponding to finite resolution results, and the upper bound to results extrapolated to infinite resolution). For model HS-nod-SN-high-accr, we only report finite resolution results (i.e. a lower bound) as the extrapolation was not provided in reference Barausse et al. (2023).

Fig. 8. Prediction of the yearly rate of mergers as a function of the total source-frame mass of the binary. Thick lines show the median and shaded areas the 90% confidence region after fitting for the LF and the GW background of EPTA. The green curves correspond to the intrinsic merger rate, while the light green curves show the rate of events detectable by LISA, assuming an S/N threshold of 10. Predictions from other models are included for comparison.

Fig. 9. Prediction of the yearly merger rate as a function of the redshift (see also Fig. 8).

In Figs. 7, 8 and 9 we also compare our results for the intrinsic rate with those of other models from the literature, selecting ones that have shown good agreement with EPTA (Antoniadis et al. 2024a): models Q3-nod, Q3-d and PopIII-d from Klein et al. (2016), model HS-nod-SN-high-accr from Barausse et al. (2023), and model Loud from Izquierdo-Villalba et al. (2024). We refer to the upcoming paper by members of the LISA astrophysics working group for a more exhaustive comparison of LISA predictions in different SAMs and cosmological simulations (LISA Astrophysics Working Group, in prep.).

The models Q3-nod, Q3-d and HS-nod-SN-high-accr assume that MBHs form from heavy seeds (from the collapse of protogalactic discs Volonteri et al. 2008), while PopIII-d assumes that they form from the remnant of Pop III stars (Madau & Rees 2001). Models HS-nod-SN-high-accr and Q3-nod also feature no delays between galaxy mergers and MBH mergers (although they include delays between halo and galaxy mergers), while PopIII-d and Q3-d do include delays due to stellar hardening, gas-driven migration and triple black-hole interactions that affect MBH binaries. Moreover, in model HS-nod-SN-high-accr the influx of gas into the nuclear regions hosting MBHs is boosted by a factor ∼4, to achieve better agreement with the PTA results (Antoniadis et al. 2024a; Barausse et al. 2023).

The Loud model by Izquierdo-Villalba et al. (2024) uses the L-Galaxies SAM with a mix of heavy and light seeds to track MBH binary evolution after galaxy mergers. It considers three phases: dynamical friction, hardening, and GW emission. The dynamical friction phase of the satellite MBHs lasts for a time given by Binney & Tremaine (2008) and starts at the radius at which the tidal forces stripped 80% of the galaxy stellar mass. At the end of this phase, the MBHs form a gravitational bound binary, entering the hardening and GW inspiral phases whose separation and eccentricity are evolved consistently based on whether the environment is gas-rich or gas-poor. The model also accounts for triple interactions and, during the MBH binary growth, an anticorrelation between the accretion rate of the primary and secondary MBH (preferential accretion). Specifically, in the Loud model, the efficiency of MBH accretion is boosted to reproduce the GW background measured by EPTA.

Our predictions agree best with the “Q3-nod” model, which assumes that MBHs form from heavy seeds and merge soon after galaxy mergers, which aligns well with the preferred values of μ_seeds and t_delay in our analysis. The HS-nod-SN-high-accr model, which makes similar hypotheses for the seeding and the delay, also shows good agreement with our results, the main difference being its large merger rate at high redshift⁸. As values of μ_seed ∼ 10² M_⊙ are also allowed by our analyis, the PopIII-d model is also mostly compatible with our results, although it is not preferred. This confirms that LISA will play a crucial role in disentangling formation scenarios.

The Loud model predicts mergers only from z ≲ 10, in tension with our results. This discrepancy may stem from the assumption in POMPOCO that the delay between galaxy and MBH mergers is independent of galaxy mass. In our model, this additional delay beyond the DF timescale is parametrised by t_delay is applied uniformly to all mergers. MBH binaries resulting from low-mass galaxy mergers are expected to experience longer formation and coalescence timescales than those in more massive systems, however. Supporting this hypothesis, the Q3-d model, which includes a post-galaxy-merger delay, yields a redshift distribution similar to the Loud model, whereas the Q3-nod model, which omits this delay, aligns more closely with our predictions. We plan to refine this aspect in future work. Note that this would not necessarily lead to a reduction in the PTA GW background, since delays following mergers of more massive galaxies–those hosting the MBHs that dominate the background could be relatively short. In this sense, our t_delay parameter may be more reflective of the delay timescales in such massive systems. In terms of total mass, the predictions of the Loud model are generally consistent with our results across most of the mass spectrum, with the exception of the high-mass end. There, it predicts significantly more mergers involving MBH binaries with m_t ∼ 10⁹ M_⊙, as a result of enhanced accretion in the model, which was introduced to fit the GW background spectrum. This is compensated in our model by a larger number of mergers occurring at low redshift, allowing us to obtain a good match with the EPTA measurement.

Regarding LISA detections, our best estimate of a few hundred events per year (see Fig. 7) is most consistent with the predictions of the Q3-nod model (∼100–200 yr⁻¹) and the PopIII-d model (∼50–100 yr⁻¹). In contrast, the predicted detection rates of ∼10 yr⁻¹ (Loud) and 10–20 yr⁻¹ (Q3-d) fall in the lower tail of our distribution, while the HS-nod-SN-high-accr model, with its prediction of ∼2000 yr⁻¹, lies in the upper tail. Detection rates are taken from Table 1 of Barausse et al. (2023) for the Q3-nod, Q3-d, PopIII-d, and HS-nod-SN-high-accr models, and from Table 1 of Izquierdo-Villalba et al. (2024) for the “Loud” model.

This comparison demonstrates that, even in its current version, POMPOCO is flexible enough to describe a range of MBH populations – a key goal of the model. To fully exploit the potential of LISA to distinguish between different formation mechanisms, however, it will be necessary to explicitly include more than one formation channel in our model in the near future.

7. Discussion and conclusions

The number of observational probes into the population of MBHs in the Universe has recently increased rapidly, and even more probes are expected in the near future. Preliminary results from the JWST already reveal intriguing discrepancies with previously accepted models, and it is anticipated that the upcoming observations will provide deeper insights. Additionally, PTA collaborations are on the verge of confirming the detection of the GW background in the nanohertz frequency band, generated by inspiralling MBH binaries. In the next decade, the launch of LISA will revolutionise our understanding of MBH formation and evolution by detecting a so far unobserved population of merging MBHs. To fully capitalise on these upcoming observations, it is essential to develop a framework that can comprehensively compare theoretical predictions with a wide variety of observational data.

To do this, we have developed the model POMPOCO, which stands for Parametrisation Of the Massive black hole POpulation for Comparison to Observations. It features 12 free parameters that are designed to effectively describe the formation and evolution of MBHs within DM halo merger trees. The computational efficiency of POMPOCO, compared to full cosmological simulations and SAMs allows us to explore a wide range of parameter space and to identify the configurations that fit the selected datasets best.

In this study, we used POMPOCO to demonstrate the consistency between the LF and the amplitude of the PTA GW background, in particular, for an enhanced LF at high redshift, as suggested by extrapolating the results from Lacy et al. (2015) and by the preliminary findings from the JWST (Harikane et al. 2023; Maiolino et al. 2024; Greene et al. 2024; Matthee et al. 2024; Taylor et al. 2025; see also Somalwar & Ravi 2025; Ellis et al. 2024b; Padmanabhan & Loeb 2024 for a similar analysis. The authors found a consistency between PTA observations and JWST black hole candidates). To achieve this, we defined a likelihood function that accommodated the LF at different redshifts between the comprehensive fits from Shen et al. (2020) and the fits to mid-IR observations from Lacy et al. (2015), and alongside the amplitude of the GW background inferred for the first frequency bin of the EPTA data. We then conducted a full Bayesian analysis. While our inferred GW spectrum clearly agrees with the EPTA data, the median of our prediction lies slightly below EPTA’s one, which might indicate a lower background, as proposed by Goncharov et al. (2024). Although we focused on the case of circular binaries in this work, we hope to include the effect of eccentricity when fitting to PTA data in future studies.

We verified our results against the observed M_* − M_BH relation at z = 0, for which we did not explicitly fit. The results agree very well with the fit by Greene et al. (2020). Furthermore, our model is capable of reproducing even some of the MBHs that were observed outside of the main distribution. Many models struggle with this feature (Habouzit et al. 2021).

A key strength of POMPOCO lies in its ability to provide a posterior distribution for the model parameters when simultaneously fitting to the LF and the GW background. For instance, our Bayesian analysis showed that approximately 10% of haloes with a mass ≳10⁷ M_⊙ at z ∼ 20 should host an MBH seed, with characteristic masses ranging from 10³ and 5 × 10⁶ M_⊙. Additionally, we estimated that a fraction of MBHs undergoes super-Eddington accretion in burst episodes lasting a few dozen million years. At cosmic noon, we predict an AGN fraction of 4 − 10% and that the peak of accretion occurs at z ∼ 2.5. Finally, we found that MBHs with masses ≳10⁸ M_⊙ experience accretion suppression by z ∼ 2.5 and that this suppression extends to MBHs with masses ≳6 × 10⁶ M_⊙ by z = 0. These estimates are consistent with the findings of more detailed numerical simulations and observations. This suggests that POMPOCO captures the essential features that drive the MBH formation and evolution.

Our model predictions for LISA are broad, but agree better with some scenarios than with others. In particular, they are more consistent with SAMs that assume that MBHs originate from heavy seeds and merge efficiently following galaxy mergers, such as the Q3-nod model (Klein et al. 2016) and HS-nod-SN-high-accr (Barausse et al. 2020). Discrepancies with the Q3-d (Klein et al. 2016) and Loud Izquierdo-Villalba et al. (2024) models in the distribution of merger redshifts suggest the need of introducing a dependence of the post-galaxy-merger delay on the mass of the host galaxy. Overall, POMPOCO is flexible enough to describe a variety of scenarios for merging MBHs while remaining compatible with current observations. To fully leverage the potential of LISA to distinguish between different formation scenarios, however, we plan to increase the flexibility of the seed distribution in our model.

Looking ahead, we aim to introduce several improvements to POMPOCO. In addition to incorporating eccentricity and other formation channels, we plan to model the spin evolution of MBHs, which is critically sensitive to the environment of merging MBH binaries (Berti & Volonteri 2008; Sesana et al. 2014; Spadaro et al. 2025). The inclusion of eccentricity and spins will also allow us to correctly model the impact of kicks following MBH mergers, which we neglected in this work, and which could lower the LISA rates. In this work, we fixed the radiative efficiency to 0.1, whereas it is expected to depend on the accretion rate and on the spin of the BHs (e.g. Merloni & Heinz 2008; Madau et al. 2014). When spins and the suppression of the radiative efficiency at very low and very high Eddington ratios are included in our model, it will properly account for this. Moreover, we intend to adopt a more physically motivated accretion model, such as one based on the mass reservoir available for accretion. This might reduce the parameter space and prevent the formation of unrealistically large MBHs. Finally, we also plan to incorporate a dependence of the time delay between galaxy mergers and the formation of MBH binaries on the masses of the host galaxies. It will be crucial to address these issues for fully extracting the astrophysical information encoded in the LISA observations and realising the immense potential of the observatory (Gair et al. 2011; Klein et al. 2016; Bonetti et al. 2019; Toubiana et al. 2021; Fang & Yang 2023; Langen et al. 2025; Spadaro et al. 2025).

Acknowledgments

We are thankful to S. Grunewald for his precious help in making POMPOCO fit to run on the Hypatia cluster, to L. Speri for providing us the EPTA data, and to M. Habouzit for fruitful discussion. A.T. is supported by ERC Starting Grant No. 945155–GWmining, Cariplo Foundation Grant No. 2021-0555, MUR PRIN Grant No. 2022-Z9X4XS, MUR Grant “Progetto Dipartimenti di Eccellenza 2023-2027” (BiCoQ), and the ICSC National Research Centre funded by NextGenerationEU. L.S. acknowledges support from the UKRI guarantee funding (project no. EP/Y023706/1). M.V. and S.B. acknowledge funding from the French National Research Agency (grant ANR-21-CE31-0026, project MBH_waves). This work has received funding from the Centre National d’Etudes Spatiales. E.B. acknowledges support from the European Union’s H2020 ERC Consolidator Grant “GRavity from Astrophysical to Microscopic Scales” (Grant No. GRAMS-815673), the European Union’s Horizon ERC Synergy Grant “Making Sense of the Unexpected in the Gravitational-Wave Sky” (Grant No. GWSky-101167314), the PRIN 2022 grant “GUVIRP – Gravity tests in the UltraViolet and InfraRed with Pulsar timing”, and the EU Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie Grant Agreement No. 101007855. This work has been supported by the Agenzia Spaziale Italiana (ASI), Project n. 2024-36-HH.0, “Attività per la fase B2/C della missione LISA”.

References

Appendix A: Halo-black hole mass relation

We show our M_h − M_BH relation at z = 0 in Fig. A.1. Unlike in the M_* − M_BH relation, we observe no break at small MBH masses, reinforcing that the break in the latter originates from the M_h − M_* relation of Moster et al. (2010).

The dotted line shows the extrapolation of the median at M_h ≲ 10¹⁴ M_⊙ to larger halo values, obtained by performing a linear fit in the range 10¹³ M_⊙ ≤ M_h ≤ 10¹⁵ M_⊙. In order to estimate the impact of underestimating MBH masses at this end on the amplitude of the GW background (given by Eq. (19)), we have have rescaled the masses of MBHs in haloes with M_h ≥ 10¹⁴ M_⊙ by the ratio between the extrapolated median and our median. We find that this increases the background by ∼5%. This results in an increase in log₁₀h_c[1/(10 yr)] of ∼0.03 dex, which is ten times smaller than our statistical uncertainty. Moreover, the difference between the median reported by EPTA and ours is ∼0.2 dex, so this potential underestimation of the MBH masses in the most massive galaxies is likely not the reason behind the discrepancy. This is because the most massive galaxies, where the difference between the extrapolated median and ours is larger, are rare and contribute little to the background. In our model, their Press-Schechter weight suppresses their contribution.

Fig. A.1. Relation between the halo mass and MBH mass at redshift z = 0.

All Tables

All Figures

Fig. 1. Evolution of the luminosity function. In green, we show the prediction of our model when fitting for the LF itself at redshifts 0.25, 0.5, 3 and 6, as well as for the GW background measured by EPTA, displaying the median (green line) and 90% confidence interval (green band). This should be compared with the fits to the observed LF from Lacy et al. (2015) (pink) and Shen et al. (2020) (purple) and the region in between (light pink), where we allow the observed LF to lie in our likelihood. At redshifts z = 4, 5 we also plot constraints from JWST (squares) from Matthee et al. (2024). The lowest luminosity point is shown in grey and not in black, because it is likely affected by incompleteness (see the discussion in Matthee et al. 2024). At redshift z = 7, we plot a lower bound obtained by Greene et al. (2024) using JWST results.

In the text

	Fig. 2. Stochastic GW background predicted by our model (green) and the free spectrum measured by EPTA (blue). The 90% confidence intervals centred at the median for the amplitude of the GW background in the first and last bins, i.e., at frequencies 1/(10 yr) and 1/(1 yr), are $4.8_{-2.7}^{+5.1}\times {10}^{-15}$ and $1.0_{-0.6}^{+1.1}\times {10}^{-15}$, respectively.
In the text

	Fig. 3. Relation of the stellar mass and MBH mass at redshift z = 0. We compare our results, in green, with the fit of the local relation from Greene et al. (2020) (‘all’), in blue, and the data also compiled in Greene et al. (2020), in black. Green points show a (random) subsample of MBHs produced in our simulations.
In the text

Fig. 4. Relation of the MBH mass and stellar mass at redshift z = 3, 4, 5, and 6. The darker green regions represent our results when restricting to MBHs with luminosity larger than 1011 L⊙, in order to mimic selection effects at high redshift. For z = 4, 5 and 6, we show JWST observations in the redshift range z = 4 to 7 (Maiolino et al. 2024; Übler et al. 2023; Harikane et al. 2023; Li et al. 2025a), some of which were used to obtain the fit of Pacucci et al. (2023) (light orange, excluding the more recent data of Li et al. 2025a). For comparison, we also show the fit to the observed local relation (light blue) of Greene et al. (2020).

In the text

	Fig. 5. Posterior on the model parameters related to seeding and BH mergers when fitting for the LF and GW background. We show the mean and standard deviation of the log-normal distribution of seed BH masses, μseed and σseed; the minimum mass of haloes seeded Mh, seed and the seeding probability fseed; the delay of binary BH mergers (in addition to halo dynamical friction), tdelay.
In the text

Fig. 6. Posterior on the model parameters related to accretion when fitting for the LF and the GW background. We show the slope of the power-law distribution γburst and the duration tburst of burst accretion; the mean and standard deviation of the log-normal distribution of steady accretion rates, μsteady and σsteady; the redshift zcut below which accretion is shut off for heavier MBHs; the slope of the increase of the shut-off mass, mcut, 0 at z = 0, with redshift, αcut.

In the text

Fig. 7. Prediction of the yearly rate of MBH mergers after fitting for the LF and the GW background measured by EPTA. The solid curve shows the intrinsic rate, the dotted one the rate of detetable events, assuming an S/N threshold of 10. For comparison, we also report the predictions for the intrinsic rate of other models from the literature (see text for description). For Q3-nod, Q3-d and PopIII-d, the shaded areas show the range reported in reference Barausse et al. (2023) (the lower bound corresponding to finite resolution results, and the upper bound to results extrapolated to infinite resolution). For model HS-nod-SN-high-accr, we only report finite resolution results (i.e. a lower bound) as the extrapolation was not provided in reference Barausse et al. (2023).

In the text

Fig. 8. Prediction of the yearly rate of mergers as a function of the total source-frame mass of the binary. Thick lines show the median and shaded areas the 90% confidence region after fitting for the LF and the GW background of EPTA. The green curves correspond to the intrinsic merger rate, while the light green curves show the rate of events detectable by LISA, assuming an S/N threshold of 10. Predictions from other models are included for comparison.

In the text

	Fig. 9. Prediction of the yearly merger rate as a function of the redshift (see also Fig. 8).
In the text

	Fig. A.1. Relation between the halo mass and MBH mass at redshift z = 0.
In the text