© The Authors 2024

1. Introduction

The large fraction of massive galaxies hosting massive black holes (MBHs) together with the hierarchical assembly of galaxies suggest that massive black hole binaries (MBHBs) are unavoidable products of galaxy evolution. The formation of these systems and their coalescence within the Hubble time is regulated by several processes (Begelman et al. 1980). Large-scale (∼kiloparsecs) dynamical friction, acting after galaxy mergers, causes the sinking of the MBHs towards the centre of the newly formed galaxy (Milosavljević & Merritt 2001; Yu 2002; Mayer et al. 2007; Callegari et al. 2009, 2011; Bortolas et al. 2022; Li et al. 2022; Amaro-Seoane et al. 2023). At smaller separations (< parsecs), the interaction with single stars and a putative circumbinary disc surrounding the two MBHs brings the two objects down to a distance where the emission of gravitational waves (GWs) lead to their final coalescence (Quinlan & Hernquist 1997; Sesana et al. 2006; Vasiliev et al. 2014; Sesana & Khan 2015; Escala et al. 2004, 2005; Dotti et al. 2007; Cuadra et al. 2009; Biava et al. 2019; Bonetti et al. 2020; Franchini et al. 2021, 2022).

Gravitational waves (GWs) emitted by MBHBs during the late inspiral, merger, and ringdown phase at the end of their lives are the main target of current and future experiments. On the one hand, the Laser Interferometer Space Antenna (LISA, Colpi et al. 2024), planned to be launched in 2035, will detect GWs at 0.1 − 100 mHz, probing MBHBs in the 10⁴ − 10⁷ M_⊙ range. On the other hand, Pulsar Timing Array (PTA) experiments seek the nHz GWs produced by cosmologically nearby (z < 1) and slowly inspiraling more massive (> 10⁸ M_⊙) MBHBs. Currently, five main PTA experiments are taking data: the European Pulsar Timing Array (EPTA, Kramer & Champion 2013; Desvignes et al. 2016), the North American Nanohertz Observatory for Gravitational Waves (NANOGrav, McLaughlin 2013; Arzoumanian et al. 2015), Parkes Pulsar Timing Array (PPTA, Manchester et al. 2013; Reardon et al. 2016), the Indian Pulsar Timing Array (InPTA, Joshi et al. 2022) and Chinese Pulsar Timing Array (CPTA, Lee 2016) projects. The recent results reported by these collaborations provide evidence of a stochastic GW background (sGWB) with amplitude A that ranges between [1.7 − 3.2]×10⁻¹⁵, at a reference frequency f  =  1 yr⁻¹ (Agazie et al. 2023; EPTA Collaboration 2023; Reardon et al. 2023; Xu et al. 2023). Different theoretical studies have been carried out to interpret the nature of the sGWB. Despite the signal being compatible with a population of MBHBs (EPTA Collaboration & InPTA Collaboration 2024; Agazie et al. 2023), modern, sophisticated hydrodynamical simulations and semi-analytical models, aiming at reproducing a wide array of cosmological observations, tend to produce smaller sGWBs, with typical amplitudes of A ≈ 1 × 10⁻¹⁵ (see e.g. Sesana et al. 2009; Kelley et al. 2017a,b; Bonetti et al. 2018a; Izquierdo-Villalba et al. 2022; Curyło & Bulik 2024; Saeedzadeh et al. 2024; Li et al. 2024). Specifically, EPTA Collaboration & InPTA Collaboration (2024) showed that state-of-the-art semi-analytical models require important changes to reach the current sGWB reported by PTA collaborations. For instance, a fast dynamical evolution for MBHBs, a rapid and larger mass growth of MBHs, and a bigger normalization in the scaling relations are fundamental requirements in semi-analytical models to reconcile theoretical predictions and current observational constraints. However, these requisites are not easy to reach unless breaking current observational constraints of the MBH population, such as the black hole mass function or quasar luminosity function (see e.g. Izquierdo-Villalba et al. 2022; Sato-Polito et al. 2023).

Due to its high amplitude, the signal reported by the PTA collaborations is already providing valuable information about the possible formation and evolution of > 10⁸ M_⊙ MBHBs. In this paper, we investigate whether this signal has also implications for LISA. In particular, we study how the forecasts for LISA MBHBs with a potential observable electromagnetic counterpart are affected when a galaxy formation model is tuned to reproduce both the latest results on the nHz sGWB and the electromagnetic emission of MBHs. To this end, we make use of the L-Galaxies semi-analytical model which is a unique framework that includes at the same level of detail the physics involved in the assembly of galaxies, MBHs, and MBHBs (Henriques et al. 2015, 2020; Yates et al. 2021; Izquierdo-Villalba et al. 2020, 2022; Spinoso et al. 2023). Under the assumption that the whole PTA signal is coming from a low-z population of MBHBs, the merger trees of the large Millennium simulations (Springel 2005) are used to explore which conditions and model modifications are required in order to match the theoretical predictions about MBHs and MBHBs with current PTA and quasars/AGN electromagnetic constraints. The effect of the new modeling in the LISA MBHBs is explored by using the Millennium-II merger trees (Boylan-Kolchin et al. 2006) whose resolution is adequate to model smaller galaxies and MBHs down to ≈10^3 − 4 M_⊙.

The paper is organized as follows: in Sect. 2 we summarise the main physics included in L-Galaxies to trace the formation and evolution of galaxies, MBHs, and MBHBs. Furthermore, we show the difficulties of current galaxy formation models to reach the high level of nHz sGWB without mismatching the number density of active MBHs. In Sect. 3 we introduce a model able to generate an sGWB compatible with the latest PTA results, and to reduce the tension seen in the quasar bolometric luminosity function. In Sect. 4 we explore the effect of the nHz signal on the population of LISA MBHBs. In Sect. 5 we underline several caveats of the model. Finally, in Sect. 6 we summarise the main results of the paper. A Λ cold dark matter (ΛCDM) cosmology with parameters Ω_m = 0.315, Ω_Λ = 0.685, Ω_b = 0.045, σ₈ = 0.9 and h = H₀/100 = 67.3/100 km s⁻¹ Mpc⁻¹ is adopted throughout the paper (Planck Collaboration XVI 2014).

2. L-Galaxies semi-analytical model

In this section, we present the semi-analytical model named L-Galaxies (SAM, Guo et al. 2011; Henriques et al. 2015, 2020; Yates et al. 2021). In brief, L-Galaxies is a code that tracks the cosmological assembly of galaxies through a set of analytical equations solved along the assembly history of dark matter halos, as given by their respective merger tree. On top of this, the latest modifications presented in Izquierdo-Villalba et al. (2020, 2022), and Spinoso et al. (2023) enable L-Galaxies to trace the formation and evolution of single and binary MBHs. We stress that all the physics described here, together with the default values of the free parameters, constitutes the SAM fiducial model tagged as Quiet model (see Table 1).

2.1. Dark matter

L-Galaxies is a flexible semi-analytical model run on top of different dark matter (DM) merger trees extracted from N-body DM-only simulations. In particular, its performance has been tested in the Millennium and TNG-DARK suit of simulations (see e.g. Henriques et al. 2015; Ayromlou et al. 2021). In this work, we use the merger trees extracted from the Millennium (MS, Springel 2005) and Millennium-II (MSII, Boylan-Kolchin et al. 2009) simulations whose minimum particle mass and large cosmological volumes allow to trace the assembly of halos in a broad mass range (10⁸ − 10¹⁴ M_⊙). MS follows the cosmological evolution of 2160³ DM particles of mass 8.6 × 10⁸ M_⊙ h⁻¹ inside a periodic box of 500 Mpc h⁻¹ on a side, from z = 127 to the present. MSII can be thought of as a high-resolution version of the MS, as it follows the same number of particles with a mass resolution 125 times higher (6.885 × 10⁶ M_⊙ h⁻¹) in a box 125 times smaller (100 Mpc h⁻¹). MS and MSII were stored at 63 and 68 epochs or snapshots, respectively. All the structures formed in these simulations were found by applying friend-of-fiend and SUBFIND algorithms and arranged with the L-HALOTREE code in the so-called merger trees (Springel et al. 2001). Given the coarse time resolution offered by the outputs of MS/MSII (snapshots are separated by ≈300 Myr), L-Galaxies performs an internal time interpolation of 5−50 Myr (depending on redshift) to improve the tracing of the baryonic physics involved in galaxy evolution. Finally, both simulations, originally run with the WMAP1 and 2dFGRS concordance cosmology, were re-scaled with the procedure of Angulo & White (2010) to match the cosmological parameters provided by Planck first-year data (Planck Collaboration XVI 2014).

2.2. Gas and stars

L-Galaxies includes a sophisticated galaxy formation model, which is able to track the evolution of the gas and stellar components of structures forming along the DM merger trees. As soon as a DM halo collapses, the model assigns to it an amount of baryons consistent with the cosmological baryonic fraction (White & Frenk 1991). These baryons are initially distributed in a quasi-static hot gas atmosphere which is able to cool down at a rate that depends on the redshift and mass of the hosting DM halo (Guo et al. 2011). The gas that is cooled falls at the centre of the DM halo leading to the formation of a gas disc due to angular momentum conservation. Star formation events occur as soon as the gas disc exceeds a critical mass, giving rise to a stellar component distributed in a disc (Croton et al. 2006). As a consequence of star formation, massive and short-lived stars explode as supernovae injecting energy and metals into the cold gas disc, reheating it, and eventually pushing it beyond the virial radius of the DM halo (Guo et al. 2011). The cold gas component of massive galaxies is also regulated by the radio-mode feedback of the central MBH, which efficiently reduces or even suppresses cooling flows (Croton 2006). Galaxies can also trigger star formation events through interaction with companion satellites. These events are divided into major and minor mergers. The first ones take place between galaxies with baryonic masses differing by less than a factor of 2 and the final result is the transformation of the remnant galaxy into a pure bulge. On the other hand, minor mergers occur during more extreme mass ratios and they are able to trigger the formation of bulges without destroying the stellar disc of the most massive galaxy. In addition to these two merger treatments, the model includes the prescription of smooth accretion presented by Izquierdo-Villalba et al. (2019) to deal with the physics of extreme-minor mergers. During these events it is expected that the stellar remnant of the satellite gets diluted inside the disc of the central galaxy before being able to reach the nucleus, thus, losing the possibility of making the bulge of the primary galaxy grow. Indeed, Izquierdo-Villalba et al. (2019) showed that the inclusion of these processes is important to recover the observed morphology of dwarf galaxies (M_stellar ≤ 10⁹ M_⊙) in L-Galaxies when this one is run on top of the MSII simulation. Besides mergers, massive disc-dominated galaxies are allowed to develop galactic bulges through disc instability (Efstathiou et al. 1982). Finally, L-Galaxies models large-scale effects such as ram pressure stripping or galaxy tidal disruption (Henriques & Thomas 2010; Guo et al. 2011).

2.3. Massive black holes

In this section, we briefly describe the main physics included in L-Galaxies to deal with MBHs.

2.3.1. Formation of massive black holes

L-Galaxies includes a refined physical model to follow the genesis of MBHs (Spinoso et al. 2023). Specifically, it tracks the spatial variations of metals and Lyman Werner radiation to account for the formation of massive seeds (10³ − 10⁵ M_⊙) via the direct collapse of pristine massive gas clouds and the collapse of dense, nuclear stellar clusters originated by early star-formation episodes. On the other hand, the formation of light seeds (10 − 100 M_⊙) after the explosion of the first generation of stars (also known as PopIII stars) is also accounted for by using a subgrid approach which takes as an input the results of GQD Press–Schechter based SAM (Valiante et al. 2021). Notice that this seeding model is only used in the MSII trees since their high halo resolution allows us to track self-consistently the genesis of the first MBHs. On the contrary, the mass resolution of MS merger trees hinders the possibility of using the seeding model presented in Spinoso et al. (2023). Therefore, when employing the MS simulation, we assume that all the newly resolved halos (which have mass of ∼10¹⁰ M_⊙), regardless of redshift, host central MBHs with a fixed mass of 10⁴ M_⊙. We have checked that the specific value of the initial MBH mass in the MS merger trees has a marginal effect in the sGWB at nHz frequencies. The cause is the fact that PTA MBHBs are placed in massive galaxies that undergo an intense merger history, feeding MBHs with large amounts of gas that erase any memory of the initial seed mass in a few hundred Myr.

2.3.2. Growth of massive black holes in the Quiet model

As soon as the MBH forms, different processes can trigger its growth. In particular, in L-Galaxies the MBH growth is divided into three different channels: cold gas accretion, hot gas accretion, and mergers with other MBHs. Among these, the first one is the main driver of black hole growth at any redshift and is triggered by galaxy mergers and disc instability (DI) events. On the one hand, after a galaxy merger, a fraction of cold gas is accreted by the nuclear black hole:

$$\begin{array}{c}\hfill \mathrm{\Delta }{M}_{\mathrm{BH}}^{\mathrm{gas}}=f_{\mathrm{BH}}^{\mathrm{merger}}{(1+z_{\mathrm{merger}})}^{\alpha }\frac{m_{\mathrm{R}}}{1+{(V_{\mathrm{BH}}/V_{200})}^2}{M}_{\mathrm{gas}},\end{array}$$(1)

where m_R ≤ 1 is the baryonic ratio of the two interacting galaxies, V₂₀₀ the virial velocity of the host DM subhalo, z_merger the redshift of the galaxy merger, M_gas the cold gas mass of the galaxy. $f_{\mathrm{BH}}^{\mathrm{merger}}$, V_BH and α are three adjustable parameters set to 280 km s⁻¹, 0.02 and 5/2, respectively, in what we refer hereafter as Quiet model (see Table 1). While $f_{\mathrm{BH}}^{\mathrm{merger}}$ and V_BH characterize the efficiency of mergers in making gas lose angular momentum and flow towards the galactic nucleus, α takes into account the fact that at high-z galaxies are more compact (Mo et al. 1998; Shen et al. 2003; van der Wel et al. 2014; Lange et al. 2015) and thus, any high-z merger event should be more efficient in bringing gas onto the MBHs compared to low-z events (see Bonoli et al. 2009; Izquierdo-Villalba et al. 2020, 2022, for similar approaches).

On the other hand, during a disc instability, the black hole accretes an amount of cold gas proportional to the mass of stars that trigger the stellar disc instability (see Izquierdo-Villalba et al. 2020):

$$\begin{array}{c}\hfill \mathrm{\Delta }{M}_{\mathrm{BH}}^{\mathrm{gas}}=max[f_{\mathrm{BH}}^{\mathrm{DI}}{(1+z_{\mathrm{DI}})}^{\alpha }\frac{\mathrm{\Delta }{M}_{\mathrm{stars}}^{\mathrm{DI}}}{1+{(V_{\mathrm{BH}}/V_{200})}^2},M_{\mathrm{gas}}],\end{array}$$(2)

where z_DI is the redshift at which the disc instability takes place and $f_{\mathrm{BH}}^{\mathrm{DI}}$ is a free parameter that takes into account the gas accretion efficiency, set to 0.0015. V_BH and α have the same value as in the case of mergers. As described in Eq. (1), all the free parameters involved in Eq. (2) try to capture the efficiency of disc instabilities in feeding with cold gas the nuclear parts of the galaxy.

After a galaxy merger or a disc instability has occurred at a time t₀, the cold gas available for accretion ($\mathrm{\Delta }{M}_{\mathrm{BH}}^{\mathrm{gas}}$) is assumed to settle in a reservoir around the black hole, M_Res¹. Instead of instantaneous gas consumption, the model considers that the gas reservoir is progressively consumed through an Eddington-limited growth phase, followed by a second phase of low accretion rates (Hopkins et al. 2005, 2006; Marulli et al. 2006; Bonoli et al. 2009). To characterize these two stages we introduce in L-Galaxies the parameter f_Edd², defined as the ratio between the bolometric (L_bol) and the Eddington luminosity (L_Edd):

$$\begin{array}{c}\hfill f_{\mathrm{Edd}}(t)=\{\begin{array}{cc}1& M_{\mathrm{BH}}(t)\le M_{\mathrm{E}}\\ \\ \frac{1}{{[1+{((t-t_0)/t_Q)}^{1/2}]}^{2/\beta }}& M_{\mathrm{BH}}(t)>M_{\mathrm{E}}.\end{array}\end{array}$$(3)

The duration of the Eddington limited phase is determined by the value of M_E = M_BH(t₀)+ℱ_EddM_Res(t₀), which is the mass reached by the MBH after consuming a fraction ℱ_Edd of its gas reservoir. Following Izquierdo-Villalba et al. (2020) the value of ℱ_Edd is set to 0.7. Note that if a galaxy undergoes a new merger or DI while the central MBH is still accreting mass from a previous event, the new cold gas driven around the MBH environment is added to the previous remnant M_Res and the growth re-starts under the new initial conditions. Finally, t_Q gives the time-scale at which f_Edd decreases and is defined as t_Q = t_d ξ^β/(βln10), with t_d = 1.26 × 10⁸ yr, β = 0.4 and ξ = 0.3. The choice of these values is based on Hopkins & Hernquist (2009) who showed that models of self-regulated MBH growth require 0.3 < β < 0.8 and 0.2 < ξ < 0.4. We highlight that any change in the values of β and ξ in the interval suggested by Hopkins & Hernquist (2009) has a small effect on our results since the bulk of the MBH growth happens during the Eddington-limited phase. Finally, during any of the events that make the MBH grow, L-Galaxies tracks the evolution of the black hole spin (a) in a self-consistent way. During gas accretion events, the model uses the approach presented in Dotti et al. (2013) and Sesana et al. (2014) which links the number of accretion events that spin-up or spin-down the MBH with the degree of coherent motion in the bulge. On the other hand, after an MBHB coalescence the final spin is determined by the expression of Barausse & Rezzolla (2009) where a distinction between wet and dry mergers is done to compute the alignment/anti-alignment between the two MBHs. For further details on the implementation of the spin model inside L-Galaxies, we refer the reader to Izquierdo-Villalba et al. (2020).

In Fig. 1 we compare the observational constraints of the quasar LFs (Hopkins et al. 2007; Aird et al. 2015; Shen et al. 2020) with the predictions of the Quiet model of L-Galaxies. As we can see, the model is in good agreement with the observations, being able to reproduce the redshift evolution of active MBHs from z = 0.5 up to z = 3.

Fig. 1. Quasar luminosity functions at z = 0.5, 1.0, 2.0, 3.0 for the Quiet (blue) and Boosted model (red) run in the Millennium merger trees (i.e. L-Galaxies + MS). The error bars correspond to the Poissonian error. The results are compared with the observations of Hopkins et al. (2007) (triangles), Aird et al. (2015) (squares) and Shen et al. (2020) (circles). The lower panel represents the ratio between the model and the data points.

2.4. Massive black hole binaries

The dynamical evolution of MBHB in L-Galaxies is divided into different stages (Begelman et al. 1980). The first one is called pairing phase and it starts after the merger of the two galaxies (see Sect. 2.2). During this process, the MBH hosted by the less massive galaxy undergoes dynamical friction which causes its sinking toward the galactic centre of its new host. The process is modeled according to the standard Binney & Tremaine (1987) equation which depends on the mass and orbital circularity of the MBH, the initial position at which the satellite galaxy deposited the MBH after the merger (∼kpc separation), and the velocity dispersion of the remnant galaxy. Notice that the model assumes that the nuclear MBH of the central galaxy is not displaced from the galactic centre after the galactic merger. This is a fair assumption since the dynamical friction time scale is dominated by the lighter MBH making irrelevant the displacement of the central and most massive MBH. Finally, we stress that the secondary MBH might be embedded inside nuclear stellar clusters (NSCs, not accounted in the model yet) leading to larger effective masses during the dynamical fiction than the ones assumed in the MBHB model of L-Galaxies. This is a caveat that will be addressed in future works by using a phenomenological model (see e.g. Polkas et al. 2023) or by including a self-consisitent NSC modelling inside L-Galaxies (Hoyer et al., in prep.).

Once the dynamical friction phase ends, the satellite MBH reaches the galactic nucleus of the new galaxy and it binds with the central MBH (∼pc separation) starting the so-called hardening phase and giving rise to a massive black hole binary. While the most massive MBH of the binary is flagged as primary, the lighter one is tagged as secondary. The initial eccentricity of the binary orbit is randomly drawn in the range [0, 0.99], while the initial semi-major axis is set to the scale at which the stellar content of the galaxy (distributed according to a Sérsic model) equals the mass of the secondary MBH. The eccentricity and separation of the MHBH are then evolved self consistently according to the environment in which the system is embedded. In case the gas reservoir around the binary (M_Res) is larger than its total mass (M_Bin) the system evolves by interacting with a circumbinary gaseous disc, following the prescription of Dotti et al. (2015). Otherwise, the system is driven by the interaction with stars according to the theory developed by Quinlan & Hernquist (1997), Sesana & Khan (2015). In both cases, GW emission takes over at smaller separations (≲mpc), leading the binary to final coalescence. We refer to Izquierdo-Villalba et al. (2022) for a detailed description of the equations used to evolve the MBHB eccentricity and separation. In the case of repeated mergers, a third MBH can reach the nucleus of the remnant before the pre-existent binary completes its evolution. In this case, an MBH triplet forms and the outcome of the triple interaction is modeled according to the tabulated values of Bonetti et al. (2018a).

On top of the dynamical evolution of MBHBs, L-Galaxies allows MBHs in the pairing and hardening phase to increase their masses. Recent hydrodynamical simulations of merging galaxies with central MBHs have shown that the secondary galaxy suffers large perturbations during the pericenter passages around the central one (Callegari et al. 2009, 2011; Capelo et al. 2015; Gabor et al. 2016). Under these circumstances, the black hole of the secondary galaxy experiences accretion enhancements mainly correlated with the galaxy mass ratio. To include these findings, L-Galaxies assumes that right before the galaxy merger, the black hole of the secondary galaxy increases its gas reservoir according to Eq. (1). Once the satellite MBH is deposited in the new galaxy and starts its pairing phase, the gas reservoir is consumed according to the two-phase model described in Sect. 2.3.2. The accretion process onto the pairing MBH lasts until it consumes the total gas reservoir stored before the merger. On the other hand, gas accretion onto MBHBs has been extensively explored by different theoretical studies (D’Orazio et al. 2013; Farris et al. 2014; Moody et al. 2019; Muñoz et al. 2019; D’Orazio & Duffell 2021). The findings of these have shown that irrespective of the mass ratio of the binaries, the gas accretion onto the secondary MBH is sufficient to modify the final mass ratio of the binary, moving the initial values toward larger ones (see e.g. Farris et al. 2014; Duffell et al. 2020). Based on this picture, L-Galaxies assumes that an MBHB in the hardening phase featuring a gas reservoir progressively consumes it according to the results of Duffell et al. (2020). In brief, the accretion rate of a primary black hole (${\dot{M}}_{{\mathrm{BH}}_1}$) is fully determined by the binary mass ratio (q) and the accretion rate of the secondary black hole (${\dot{M}}_{{\mathrm{BH}}_2}$) by:

$$\begin{array}{c}\hfill {\dot{M}}_{{\mathrm{BH}}_1}={\dot{M}}_{{\mathrm{BH}}_2}(0.1+0.9q).\end{array}$$(4)

Except in the case of equal mass systems, the secondary MBHs are farther away from the binary centre of mass than primary ones. This causes them to be closer to the circumbinary disc edges and thus display high accretion rates. Based on this, L-Galaxies fix the accretion of the secondary black hole at the Eddington limit and determine the accretion onto the primary according to Eq. (4).

2.5. Model constraints from Pulsar Timing Arrays

Recent PTA results suggest the existence of an sGWB at nHz frequencies, compatible with a population of merging low-z MBHBs (Agazie et al. 2023; EPTA Collaboration 2023; Reardon et al. 2023; Xu et al. 2023)³. Despite the current significance is still below the canonical 5σ threshold, this signal provides galaxy formation models with a new condition to calibrate their underlying MBH growth physics, adding to the standard electromagnetic constraints coming from observations of galaxy properties and the quasar luminosity function.

The foundation to perform a comparison between recent PTA results and galaxy formation models resides in determining the comoving number density of MBHB mergers (d²n/dzdℳ) per unit redshift, z, and rest-frame chirp mass, ℳ⁴. Following Sesana et al. (2008) and making the specific assumption that inspiralling MBHBs in the PTA band are in circular orbits evolving purely due to GW emission, the characteristic sGWB can be written as:

$$\begin{array}{c}\hfill h_{\mathrm{c}}^2(f)=\frac{4G^{5/3}{f}^{-4/3}}{3c^2{\pi }^{1/3}}{\displaystyle \int \int }\mathrm{d}z\mathrm{d}\mathcal{M}\frac{d^{2}n}{\mathrm{d}z\mathrm{d}\mathcal{M}}\frac{{\mathcal{M}}^{5/3}}{{(1+z)}^{1/3}},\end{array}$$(5)

where f is the frequency of the GWs in the observer frame. This expression is often simplified as:

$$\begin{array}{c}\hfill h_{\mathrm{c}}(f)=A{\left(\frac{f}{f_0}\right)}^{-2/3},\end{array}$$(6)

where A is the amplitude of the signal at the reference frequency f₀. Hereafter, we will set f₀ = 1 yr⁻¹ and refer A to that frequency. Under these assumptions, the Quiet model of L-Galaxies predicts an d²n/dzdℳ which generates a sGWB of amplitude A ∼ 1.2 × 10⁻¹⁵ (see Table 1), fully ruled by z < 1 MBHBs with ℳ > 10⁸ M_⊙ (see Izquierdo-Villalba et al. 2023). Despite agreeing with past theoretical studies (Wyithe & Loeb 2003; Sesana et al. 2008, 2009, 2016; Sesana 2013; Ravi et al. 2015; Bonetti et al. 2018b; Kelley et al. 2017a; Siwek et al. 2020) its value is lower than the one reported by PTA collaborations which span over the range [1.7 − 3.2]×10⁻¹⁵ (EPTA Collaboration 2023; Agazie et al. 2023; Reardon et al. 2023; Xu et al. 2023). Interestingly, such tension is also present in other recent SAMs with MBHBs such as Li et al. (2024) or Curyło & Bulik (2024).

To reconcile theoretical predictions with the recent PTA results, EPTA Collaboration & InPTA Collaboration (2024) explored changes in the dynamical models of MBHBs included in L-Galaxies (e.g. faster dynamical friction phases and only stellar hardening) showing that these modifications are not enough to reduce the discrepancy. Despite this, in a recent paper, Barausse et al. (2023) showed that a Press–Schechter based SAM which includes a heavy MBH seeding scenario and assumes no delay between galaxy and MBH mergers would favor large sGWB, compatible with the latest PTA results. Another approach was reported in Izquierdo-Villalba et al. (2022) which proposed that MBHs should be more efficient in accreting cold gas after mergers and/or disc instabilities, thus becoming more massive when they enter the PTA frequency band. To increase the efficiency of MBH growth the authors explored in L-Galaxies a boosted model in which the two-phase growth model of Sect. 2.3.2 was untouched but the gas accretion efficiency, $f_{\mathrm{BH}}^{\mathrm{merger}}$ in Eq. (1), increases during mergers⁵ (hereafter Boosted model, see Table 1). The results showed that the Boosted model was able to generate a stochastic GW background of amplitude [1.9 − 2.6]×10⁻¹⁵, compatible with the latest PTA measurements. Despite this better agreement, the increase of the MBH growth was hindering the possibility of reproducing the number density of active MBHs. This is presented in Fig. 1 where we can see that the Quiet model is able to follow the observational trends presented in Hopkins et al. (2007), Aird et al. (2015), and Shen et al. (2020). However, the Boosted one displays a systematic overprediction (up to ∼2 dex) in the number density of bright quasars (L_bol ≳ 10⁴⁶ erg s⁻¹) at z < 2.

The results presented above suggest that increasing the mass of MBHs is a good avenue to reconcile theoretical models and PTA observations. However, the physical mechanism and the time scale by which the MBHs reach the necessary mass to generate a loud sGWB should be addressed carefully since the over-prediction of the quasar luminosity function seems to be a natural consequence. Taking this into account, in the next section we explore the effect of allowing super-Eddington accretion episodes to grow the population of single MBHs. Specifically, we investigate if these events enable the assembly of a population of MBHB compatible with the nHz sGWB and generate a population of active MBHs in agreement with the quasar luminosity functions. We stress that super-Eddington accretion is only enforced on single MBHs while the growth of MBHBs is modeled in the same way as described in Sect. 2.4.

3. A faster assembly for the MBH population

In this section, we explore the possibility of extending the Quiet model of L-Galaxies in such a way that galaxy mergers are more efficient in fuelling gas onto MBHs and some MBHs, under certain conditions, can undergo super-critical accretion events. While the former requirement is done in the same way as in the Boosted model (see Sect. 2.5), the latter is described in the next section. The interplay between these two processes is calibrated by running L-Galaxies on the Millennium merger trees. Finally, we stress again that super-Eddington accretion events will be allowed only to nuclear single MBHs, being the accretion onto MBHBs modelled in the same way as Sect. 2.4.

3.1. A toy model for super-Eddington growth

Super-Eddington accretion refers to growth episodes that proceed extremely rapidly, breaking the Eddington rate and increasing the MBH mass on time scales shorter than what is allowed by the standard Eddington-limited model (Abramowicz et al. 1988). Recent simulations have shown that not all the environments in which MBHs are embedded can trigger these extreme accretion events. Only dense and dusty gas environments around the MBHs replenished by large gas inflows after galactic mergers provide the ideal conditions to trigger super-Eddington growth (see e.g. Inayoshi et al. 2016; Takeo et al. 2018; Regan et al. 2019; Toyouchi et al. 2021; Sassano et al. 2023; Massonneau et al. 2023). To account for these requirements in the L-Galaxies SAM, we assume that a super-Eddington phase is active only if: (i) at the moment of the merger/disc instability (t₀) the gas reservoir around the single MBH exceeds by a factor of ℛ^th the MBH mass:

$$\begin{array}{c}\hfill \frac{M_{\mathrm{Res}}(t_0)}{M_{\mathrm{BH}}(t_0)}>{\mathcal{R}}^{\mathrm{th}},\end{array}$$(7)

and (ii) the rate at which the gas is infalling towards the galactic centre, M_infolw, overcomes a certain threshold given by

$$\begin{array}{c}\hfill M_{\mathrm{inflow}}=\frac{\mathrm{\Delta }{M}_{\mathrm{BH}}^{\mathrm{gas}}}{t_{\mathrm{dyn}}}>M_{\mathrm{inflow}}^{\mathrm{th}},\end{array}$$(8)

where $\mathrm{\Delta }{M}_{\mathrm{BH}}^{\mathrm{gas}}$ corresponds to the cold gas that fuels the new accretion event onto the nuclear single MBH and it is determined by Eqs. (1) and (2). We set the dynamical time of the gas to inflow towards the galactic nucleus, t_dyn, as $V_{\mathrm{disc}}/R_{\mathrm{gas}}^{\mathrm{sl}}$ being V_disc the maximum circular velocity of the cold gas (assuming an exponential disc profile) and $R_{\mathrm{gas}}^{\mathrm{sl}}$ the cold gas scale length radius (see Guo et al. 2011, for the detailed model of galactic sizes included in L-Galaxies).

In case both conditions are satisfied, the MBH does not follow the Eddington-limit growth of Eq. (3) but undergoes a super-critical accretion event whose lightcurve is characterized by the following f_Edd:

$$\begin{array}{c}\hfill f_{\mathrm{Edd}}(t)=\{\begin{array}{cc}B(a)[\frac{0.985}{{\dot{M}}_{\mathrm{Edd}}/\dot{M}+C(a)}& \\ +\frac{0.015}{{\dot{M}}_{\mathrm{Edd}}/\dot{M}+D(a)}]& M_{\mathrm{BH}}(t)\le M_{\mathrm{SE}}\\ \\ \frac{1}{{[1+{((t-t_0)/t_Q)}^{1/2}]}^{2/\beta }}& M_{\mathrm{BH}}(t)>M_{\mathrm{SE}}.\end{array}\end{array}$$(9)

This two-phase-lightcurve tries to mimic the fact that even in an environment favourable for super-Eddington accretion, the AGN feedback resulting from the gas accretion makes the MBH reach a self-regulated phase within a few Myr⁶ (see e.g. Massonneau et al. 2023). The parameter M_SE is the maximum mass reached by the MBH during the super-critical accretion and it is defined as M_SE = M_BH(t₀)+ℱ_SEM_Res(t₀), being M_BH(t₀) and M_Res(t₀) the mass of the MBH and the reservoir at the moment of the (major/minor) merger and/or disc instability (occurring at t₀). ℱ_SE determines the fraction of the gas reservoir consumed by the MBH throughout the super-Eddington phase, before the large energy released during the accretion swaps the gas material around the MBH and hinders a subsequent Eddington limit phase (Lupi et al. 2016; Regan et al. 2019). For simplicity, and to reduce the large number of parameters we fix ℱ_SE = 0.1. The functions B(a), C(a) and D(a) are taken from Madau et al. (2014) and they scale with the spin of the MBH (a) as B(a) = (0.9663 − 0.9292a)^−0.5639, C(a) = (4.627 − 4.445a)^−0.5524 and D(a) = (827.3 − 718.1a)^−0.7060. Notice that we do not make any assumption about the spin value of the MBH since it is computed self-consistently in L-Galaxies after any MBHB merger and gas accretion episode (see Izquierdo-Villalba et al. 2020, for further details)⁷. Finally, Ṁ and Ṁ_Edd are the accretion rate and Eddington accretion rate onto the MBH, respectively. To determine Ṁ we extract a random number between [0 − 10⁵] distributed according to Ṁ⁻¹. This choice is motivated by recent theoretical works that show a power-law decline in the distribution of simulated populations of accreting MBHs for values above the Eddington-limit (see e.g. Fanidakis et al. 2012; Griffin et al. 2019 or Shirakata et al. 2019)⁸.

3.2. Setting up the super-Eddington conditions: The Loud model

The large number of possible galaxy merger histories provided by L-Galaxies and the Millennium merger trees enable us to explore which are the most common gas environments around MBHs and the typical amount of gas flowing toward the galactic centre after any secular or merger process. To illustrate the interplay between these quantities, Fig. 2 shows the plane M_inflow versus M_Res(t₀)/M_BH(t₀) for the three different events able to trigger the MBH growth: major merger, minor mergers and disc instabilities. These quantities have been computed under the assumption of a more efficient fuelling of gas onto MBHs (see the free parameters used in the third row of Table 1). As shown, at z > 3 galaxy interactions can trigger large gas inflows (> 10 M_⊙ yr⁻¹) that bury MBHs in large gas reservoirs of up to > 10³ times more massive than the MBH itself. This trend is the result of the fact that interacting high-z galaxies are compact, gas-rich, and host MBHs whose mass is several orders of magnitude smaller than the cold gas component. On the other hand, disc instabilities occurring at high-z do not trigger the large gas inflows shown in mergers, with rates that are typically 2 orders of magnitude smaller (0.1 M_⊙ yr⁻¹). This points out that in the high-z universe, galactic mergers between gas-rich systems are the unique systems that can sustain significant gas inflows to trigger potential super-critical accretion events onto MBHs. These results align with the findings of other works such as Pezzulli et al. (2016, 2017) and Trinca et al. (2022) who showed, by using a semi-analytical model, that small MBH seeds at high-z can increase several orders of magnitude their masses trough super-critical accretion after gas-rich galactic mergers. Furthermore, similar conclusions have been drawn from the hydrodynamical simulations of Lupi et al. (2024), which pointed out that gas-rich environments at high-z redshift can support long-lasting (tens of Myr) super-Eddington accretion phases, speeding up the growth of intermediate-mass MBHs. Figure 2 also shows that at intermediate redshifts, 1 < z < 3, the picture is similar to the one seen at higher-z. However, there is an important decrease of events with large inflows and a rise of cases where MBHs are surrounded by small gas reservoirs. This is the consequence of the fact that a large fraction of the galaxy population has already transformed its gas component into stars; as a consequence, mergers between gas-rich systems are less common than at higher redshifts. Specifically, in the case of major mergers, the peak of the M_inflow distribution is displaced down to ∼0.1 M_⊙ yr⁻¹, about 1.5 orders of magnitude smaller than the higher-z case. Finally, mergers and disc instabilities occurring at z < 1 are very inefficient in sustaining large inflows towards the galactic nucleus, with rates that rarely exceed 0.01 M_⊙ yr⁻¹. Such low inflows imply that MBHs are systematically embedded in small gas reservoirs whose masses can be two orders of magnitude smaller than the central MBH.

Fig. 2. Minfolw − MRes(t0)/MBH(t0) plane in three different redshifts bins according to L-Galaxies run on the Millennium merger trees. While the y-axis depicts how gas-rich is the environment around the MBH, the x-axis illustrates how powerful are the gas inflows towards the MBH. The blue horizontal and vertical lines correspond to the chosen ℛth, and $M_{\mathrm{inflow}}^{\mathrm{th}}$. Each row corresponds to each of the events that can trigger cold gas inflows toward the MBH: major mergers (top), minor mergers (middle) and disc instabilities (bottom). The plot has been computed by using the free parameters of the Boosted model reported in Table 1. To guide the reader, major (minor) mergers correspond to galaxy interactions which involve galaxies with baryonic masses differing by less (more) than a factor of 2.

Taking into account the trends presented above and how the LFs evolve with different thresholds in $M_{\mathrm{inflow}}^{\mathrm{th}}$ and ℛ^th (presented in Appendix A for the sake of brevity), we have chosen M_inflow = 10 M_⊙ yr⁻¹ and ℛ^th = 2 × 10³ as the best set of parameters to reproduce the evolution of the quasar population. The resulting model is tagged as Loud (see Table 1) and its LFs are presented in Fig. 3. As shown, the match between predictions and observations improves with the new model. Specifically, the Loud model matches the z > 4 LFs and predicts at z ≤ 2 up to a factor of 2 less objects for any bin of luminosity larger than 10⁴⁶ erg s⁻¹, reducing the tension seen in the Boosted case. These improvements are the result of the faster MBH population assembly which consumes most of its gas at high-z and evolves quiescently in the low-z Universe (see Appendix A). This effect can be seen in the upper and middle panel of Fig. 4 which presents the number density of MBHs with mass > 10⁸ M_⊙ in the Loud and Boosted model (see the global assembly of the MBH population of the Loud model in Appendix B). As shown, the population of > 10⁸ M_⊙ is in place much earlier in the former model than in the latter, with number densities that can be up to ∼1 dex larger. Regarding the AGN activity of such population, we can see that they are mainly active (f_Edd > 0.01) at z > 3 but they become rapidly inactive (f_Edd < 0.01) at z < 2.

Fig. 3. Evolution of the quasar luminosity function produced by the Loud model (black line), Boosted model (red line) and Quiet model (blue line) when run on the Millennium merger trees. The error bars correspond to the Poissionian error. The results are compared with the observations of Hopkins et al. (2007) (orange triangles), Aird et al. (2015) (purple squares) and Shen et al. (2020) (blue circles). The panels below each luminosity function represent the ratio of the luminosity functions predicted by Loud and Boosted model. Notice that the two lower ones do not display any line because their ratio goes beyond the limit 1.

The better agreement between the Loud model and the electromagnetic constraints does not imply a small value of the sGWB amplitude which instead has a value A = 1.8 × 10⁻¹⁵, consistent with the 90% credible interval reported by all the PTA collaborations. The reason why the sGWB is higher in the Loud model with respect to the Quiet one can be seen in the lower panel of Fig. 4, which shows the number density of MBHBs with ℳ > 10⁸ M_⊙, i.e. the systems which contribute the most to the nHz sGWB signal (see Sesana et al. 2008; Sesana 2013; Izquierdo-Villalba et al. 2022). As we can see, these MBHBs are in place much earlier in the Loud model than in the Quiet one with number densities up to 2 times larger at any redshift (especially at z < 1). Finally, the implications that our new modelling has about the population of super-Eddington sources can be found in Appendix C.

Fig. 4. Number density of active and inactive MBHs/MBHBs. Upper panel: number density of MBHs with mass > 108 M⊙ as a function of redshift, z, for Loud (black) and Boosted model (grey). Middle panel: number density of active (fEdd > 0.01) and inactive (fEdd < 0.01) MBHs with mass > 108 M⊙ as a function of redshift, z, for Loud (black) and Boosted model (grey). Lower panel: number density of MBHBs with chirp mass ℳ > 108 M⊙ as a function of redshift, z, for Loud (black) model. For reference, it has been shown the results for the Quiet model (grey). In all panels, the error bars correspond to the Poisson error. We remind the reader that the comoving volume provided by the Millennium simulation corresponds to ∼3.6 × 108 Mpc3.

4. Constraining low-redshift LISA MBHB mergers from the nHz sGWB

In this section, we study the implications of the Loud model, featuring a sGWB of amplitude 1.8 × 10⁻¹⁵ at f = 1 yr⁻¹, on the expected merger rate and electromagnetic counterpart detection of LISA MBHBs. Since we are interested in a population of MBHs up to 4 orders of magnitude smaller than the PTA one, instead of using the Millennium merger trees we will make use of the Millennium-II ones, which offer the possibility of resolving halos down to ∼10⁸ M_⊙ (see Sect. 2.1). On top of this, the high-resolution offered by Millennium-II enables us to seed MBHs in newly formed galaxies according to the multi-flavour seeding model of Spinoso et al. (2023).

4.1. The detection rate of LISA MBHBs

The redshift versus total rest-frame MBHB mass plane of merging binaries in the Loud model is presented in Fig. 5. As shown, at M_Bin < 10⁶ M_⊙ the Loud model displays many mergers occurring at z > 3 but the large majority of the events lie at 1 < z < 3 (see similar results from hydrodynamical simulations in Salcido et al. 2016). For more massive systems, most mergers occur at z < 1. To explore how the PTA signal affects the LISA predictions, in Fig. 6 we present the M_Bin − z plane where each pixel represents the ratio between the number of MBH mergers predicted by the Loud and Quiet model. While the former displays an sGWB of A = 1.8 × 10⁻¹⁵, compatible with PTA results, the latter produces a signal with A = 1.2 × 10⁻¹⁵ and is in line with a large number of past theoretical works (see e.g. Jaffe & Backer 2003; Sesana et al. 2008; Sesana 2013; Roebber et al. 2016; Kelley et al. 2017a; Barausse et al. 2020). As shown, the Loud model produces 2 − 5 times more mergers of MBHBs with total masses > 10⁶ M_⊙ than the Quiet one. These differences decrease for 10⁵ < M_Bin < 10⁶ M_⊙ where the ratio varies between 0.9 − 1.5. For the lighter systems (< 10⁴ M_⊙) the Loud model displays a deficit with respect to the Quiet one. Specifically, the number of MBHB mergers at these masses can be down to 0.5 less frequent at z < 4. The decrease of light MBHB mergers in favour of more massive ones is the result of an earlier and faster assembly of the MBH population in the Loud model than in the Quiet one (see Fig. 4). The different MBHB populations generated by the two models point towards different merger rates. These are summarised in Table 2. As shown, the Loud and Quiet model displays a similar global rate of 12.7 yr⁻¹. However, as expected from Fig. 6 the differences are more evident when dividing the population by masses. Specifically, the Loud model predicts merger rates than are ∼1.2 smaller than the Quiet case when MBHBs of total mass 10^3 − 5 M_⊙ are considered. However, it boosts by a factor of 1.5 the coalescences involving systems with total mass > 10⁵ M_⊙ (0.8 − 4.3 yr⁻¹, depending on the mass).

Fig. 5. Distribution of MBHB mergers predicted by the model Loud in the plane redshift (z) versus total rest-frame MBHB mass (MBin). Each pixel is encoded by the median signal-to-noise ratio (S/N). The contours represent the number of objects in the z − MBin plane. The small plots around the z − MBin plane correspond to the redshift (right panel) and source-frame mass (top panel) distribution of the differential number of MBHBs. In the right panel, the total distribution has been divided into 5 different mass bins.

Having determined the MBHB merger rate, it is important to calculate how many events are detectable by LISA. To this end, we compute the signal-to-noise ratio (S/N) of all the simulated binaries, defined as:

$$\begin{array}{c}\hfill S/N={{\displaystyle \left({\int }_{f_0}^{f_{\mathrm{f}}}{\left[\frac{h_{\mathrm{c}}(f_{\mathrm{obs}}^{\prime })}{h_{\mathrm{n}}(f_{\mathrm{obs}}^{\prime })}\right]}^2\frac{\mathrm{d}{f}_{\mathrm{obs}}^{\prime }}{f_{\mathrm{obs}}^{\prime }}\right)}}^{1/2},\end{array}$$(10)

where h_n corresponds to the characteristic strain noise parameterized as in Babak et al. (2021), while h_c represents the characteristic strain amplitude of the source defined as $h_{\mathrm{c}}(f)=4f^2{|\stackrel{\sim }{h}(f)|}^2$ being $\stackrel{\sim }{h}(f)$ the Fourier transform of the strain signal, here computed according to the phenomenological frequency-domain gravitational waveform model PhenomC (Santamaría et al. 2010). f₀ represents the starting frequency of inspiralling binaries, set for simplicity to the low-frequency cut-off limit of the LISA sensitivity curve (10⁻⁴ Hz). Instead, f_f is the maximum frequency of the signal set to 0.15 c³/G(1 + z) M_Bin⁹. The resolution frequency bin, df, used to integrate Eq. (10) is set to df = 1/T_obs where T_obs = 4 yr corresponds to the length of LISA observations. Finally, the effective spin of the MBHBs is taken from the predictions of L-Galaxies and the eccentricity is set to 0, for simplicity.

The S/N of the MBHBs generated with the Loud model are depicted in Fig. 5. Each bin of the M_Bin − z plane encodes the median S/N value of the MBHBs falling within it. As shown in Fig. 5, MBHBs of > 10⁸ M_⊙ have S/Ns that rarely surpass a value of 10. Conversely, MBHB mergers with masses 10⁵ < M_Bin < 10⁷ M_⊙ display large S/Ns, with values that can span between 100 < S/N < 10⁴. For lighter systems (< 10⁴ M_⊙) the S/N decreases and the values are systematically < 100, being the smallest for mergers of MBHBs with total mass < 5 × 10³ M_⊙ at z > 2. Taking into account the S/N distribution of the MBHBs, Table 2 presents the LISA detection rate assuming a minimum signal-to-noise ratio of 10. As we can see, the Loud model predicts a detection rate of 12.24 MBHBs per year, being the systems with 10⁴ < M_BHBs < 10⁵ M_⊙ and 10⁵ < M_Bin < 10⁶ M_⊙ the ones with largest rates (5.10 yr⁻¹ and 4.37 yr⁻¹, respectively). As expected, MBHBs with masses > 10⁷ M_⊙ are the most affected by the LISA sensitivity curve. In this mass range, LIA can detect only half of the total events (0.4 yr⁻¹). As discussed in Fig. 6, the main difference between the Quiet and Loud models resides in the fact that the latter predicts larger merger events only for MBHBs of > 10⁵ M_⊙. Since LISA is not very sensible to MBHBs of mass > 10⁷ M_⊙ it implies that the overall detected rate predicted by Quiet and Loud model does not differ much: 12.25 yr⁻¹ versus 12.24 yr⁻¹.

Fig. 6. Distribution of the z − MBin when each pixel is encoded by the ratio of the number of MBHB mergers in the Loud and Quiet model.

The results presented in this section highlight how the recent PTA signal is not significantly informative about the expected global merger rate which will be inferred from LISA observations (see opposite conclusions in Steinle et al. 2023). However, under the assumption that an MBHB population is entirely responsible for the signal observed by PTAs, our detailed galaxy formation model suggests that the loud PTA signal would favour a more numerous population of MBHs and MBHBs of > 10⁸ M_⊙ at low-z that usually expected (see Figs. 4 and B.1) which increases the possibilities (1.5 times larger than expected) of LISA to detect such kind of events during its lifetime. However, the physical processes involved in the efficient growth of large MBHs do not act in the same way for the smallest MBHs (< 10⁵ M_⊙) which in turn are the preferred targets of LISA. These objects are not able to increase their mass as fast as the most massive MBHs because the small galaxies where they are hosted are not capable of sustaining large and continuous gas inflows towards their centres after galaxy interactions. As a result, the population of low-mass MBHs (and consequently low-mass MBHBs, < 10⁵ M_⊙) is only mildly affected by the modifications of MBH growth presented in this work. Thus, small changes are seen in the global population of LISA MBHB mergers (dominated by low-mass MBHs) with respect to the Quiet model.

4.2. Electromagnetic emission of low-z LISA MBHBs

Our results suggest that the expected merger rate detected by LISA is not significantly affected by the specific level of the nHz sGWB, as shown in Table 2. Despite this, Fig. 6 shows that a model in agreement with current constraints on the nHz sGWB predicts that the LISA detection rate of MBHBs above 10⁵ M_⊙ at z < 3 could be larger compared to a model with a smaller amplitude of the sGWB. The detection and merger rates at z < 3 are depicted in the middle panel of Table 2. As shown, for MBHBs of M_Bin > 10⁵ M_⊙ Loud model displays up to a factor [1.3 − 1.5] larger rates than the Quiet one. Since these LISA MBHBs are the ones with the brightest electromagnetic emission, this implies that during the lifetime of the LISA mission, there will be higher chances to detect the electromagnetic emission of MBHBs than previously estimated. To explore the electromagnetic detectability of LISA systems in the Loud model, Fig. 7 presents the cumulative distribution function of the bolometric (and hard X-rays) luminosity of detectable MBHBs with M_Bin > 10⁶ M_⊙. The figure shows that these systems have a 50% (20%) probability of shining at L_bol > 10⁴³ erg s⁻¹ (L_bol > 10⁴⁴ erg s⁻¹). In X-rays (2 − 10 keV, see Merloni et al. 2004 for the bolometric correction)¹⁰ we see similar trends, with 50% of the MBHB with M_Bin > 10⁵ M_⊙ shining at an X-ray luminosity > 10⁴² erg s⁻¹. These prospects imply that future X-ray observatories such as Athena (Nandra et al. 2013) or optical surveys such as the Vera C. Rubin Observatory (Ivezić et al. 2019) which feature low limiting fluxes will be able to detect the electromagnetic emission coming from several LISA MBHBs.

Fig. 7. Cumulative distribution function (CDF) of bolometric (black) and hard X-ray (2 − 10 KeV, red) luminosity of LISA detectable (S/N > 10) MBHBs with masses MBin > 105 M⊙. Darker lines correspond to the Loud model whereas the light curves depict the predictions of the Quiet one. Horizontal lines highlight the CFR values of 0.5 and 0.2. The vertical pink line highlights the luminosity value of 1043 erg s−1, while black dashed and dotted lines represent the minimum luminosity that a source must have respectively at z = 1 and z = 0.5 to be detected by Athena X-ray observatory assuming a flux limit in the hard X-ray of 2 × 10−15 erg s−1 cm−2 (see Lops et al. 2023).

Having seen that around 50% of the MBHBs detected by LISA will have an observable electromagnetic counterpart, it is interesting to determine the detection rate of these multi-messenger MBHBs. To this end, the lower part of Table 2 summarises for the Loud model the detection rate of MBHBs at z < 3 featuring a bolometric luminosity > 10⁴³ erg s⁻¹. The results show that LISA will be able to detect that type of system at a rate of 2.41 yr⁻¹ (against the 2.13 yr⁻¹ displayed by the Quiet model). When the population is divided into bins of mass, we can see that MBHBs with 10⁵ − 10⁶ M_⊙ have a detection rate of 1.5 yr⁻¹ whereas MBHBs of 10⁶ − 10⁷ M_⊙ have a rate up a factor 2 smaller. For the case of MBHBs with M_bin > 10⁷ M_⊙, the Loud model shows that the detection rate can be up to 0.22 yr⁻¹ (a factor 2 larger than in the Quiet model).

5. Caveats

In this section, we discuss several caveats to take into consideration when interpreting the results presented in this work.

5.1. Multiple avenues to enlarge the nano-Hz stochastic GWB

In the recent work, EPTA Collaboration & InPTA Collaboration (2024) explored whether state-of-the-art galaxy and MBH formation and evolution models could reproduce the nHz signal reported by the EPTA collaboration. Under the assumption that the whole signal is generated by an MBHB population, the authors showed that under standard assumptions on the MBH and MBHB evolution, SAM models generally predict a signal approximately a factor of two smaller than what is detected. To overcome this limitation, it was shown that a faster dynamical evolution of MBHBs after a galaxy merger and/or a rapid and larger growth of MBHs should be required. In this work, we have explored the latter, by studying the possibility that a quick assembly of the MBH population could be caused by a larger efficiency of galactic mergers in bringing gas towards the galactic nucleus, triggering super-Eddington accretion events onto single MBHs. However, this approach is not the unique avenue to reconcile a high sGWB amplitude with galaxy formation models. For instance, the dynamical models for MBHBs, generally relying on oversimplified assumptions, could be revisited. The changes in the dynamics of MBHBs would imply variations in the population of MBHBs that could enlarge the nHz sGWB without imposing any change in the population of single MBHs. Therefore, different approaches used to reach large nHz sGWB would imply different consequences for low-z LISA MBHBs than those found in this work. In a future paper, we plan to revisit the dynamical and growth model of MBHBs, exploring which of the requirements involved in these processes increase the nHz signal.

Finally, the results presented in this work do not follow the same trends as the ones reported in Barausse et al. (2023). By using a Press–Schechter based SAM calibrated against the PTA results, the authors showed that LISA forecasts are strongly affected by the underlying PTA signal. Specifically, Barausse et al. (2023) reported that no time delays between galactic and MBH mergers, higher accretion rates onto MBHs, and heavy MBH seeding scenarios would favour large nHz sGWBs with LISA detection rates varying between 3 yr⁻¹ up to 9600 yr⁻¹. In the case of L-Galaxies, a model without delays causes a decrease in the PTA signal as the result of the smaller chirp masses of MBHBs at the time of merger (h_c ∝ ℳ^5/3, see Fig. 5 of Izquierdo-Villalba et al. 2022): the time spent by large MBHBs in the dynamical friction phase enlarges the mass of the satellite and central MBH, and the gas accretion during the hardening phase tends to make more equal mass systems, increasing in this way the chirp mass of the MBHB at merging time. The discrepancies seen between this work and Barausse et al. (2023) point out that further investigations are required to shed light on how PTA detections would impact the forecast about LISA MBHBs. Specifically, an analysis of why two SAMs that follow very accurately galaxy, MBH and MBHB evolution provide so different results would help in sharpening our knowledge about what is the main physics shaping the population of MBHBs.

5.2. LISA merger rate: Effects of the underlying assumptions

In this work, we have given predictions about the expected MBHB mergers that LISA will detect. However, several assumptions included in L-Galaxies could affect these results. One of them concerns the MBH seeding model. Whereas assumptions in the seed mass do not have any relevant impact in PTA studies, the range of masses and redshifts covered by LISA makes the results sensible to the specific physics assumed about the MBH formation. To have an idea about how the seeding model of L-Galaxies affects the LISA merger rates, we have explored a pessimistic scenario in which the formation of MBHs in L-Galaxies halos is suppressed¹¹. In this pessimistic scenario, the LISA merger rate can drop a factor of 3 (∼4 yr⁻¹) with respect to the values presented in this work.

Another caveat to take into account concerns the growth of small MBHs. Recent works pointed out that SN feedback can influence the environments around small MBHs and hinder their growth (see e.g. Habouzit et al. 2017). This scenario, not taken into account by L-Galaxies, could generate a population of lighter MBHs compared to our fiducial runs. Consequently, the dynamical friction phase could be longer and hinder the formation of MBHBs. To take into account this, we have explored a toy model in which the growth of MBHs in L-Galaxies is suppressed as long as the escape velocity of its hosting galaxy is smaller than 270 km s⁻¹ (see similar approach presented in Barausse et al. 2020). The results showed that the merger rates dropped down to 9 yr⁻¹ and the peak of the distribution was delayed a few Gyr. Further analysis on the suppression of MBH growth in dwarf galaxies will be explored in Spinoso et al. (in prep.).

Finally, the predicted LISA detection rate could be also affected by the halo mass resolution provided by the merger trees of the Millennium-II simulation. To have a tentative estimation of this, we have used L-Galaxies on top of the merger trees of the TNG50-Dark and TNG50-2-Dark simulations (Nelson et al. 2015)¹². TNG50-Dark-2 can be approximated to Millennium-II given their similar particle mass resolutions. On the other hand, TNG50-Dark can be interpreted as the high-resolution run of Millennium-II given its 10 times better mass resolution. To set an upper limit of MBHB mergers lost just by resolution we have compared the MBHB merger rates predicted by L-Galaxies on top of TNG50-Dark and TNG50-2-Dark by assuming that all the halos are seeded with an MBH of 100 M_⊙ (maximal MBH occupation) and no delays in the formation of MBHB after the galaxy merger (efficient sinking). Such comparison showed that the main differences were concentrated at z > 1 where the MBHB merger rate was a factor 1.1 larger in the TNG50-Dark than in TNG50-2-Dark (i.e. Millennium-II like run). Considering this resolution test, we can conclude that the LISA merger events discussed in this study remain almost unchanged in case the mass resolution of the underlying merger trees is raised by one order of magnitude.

6. Conclusions

In this paper, we explored if the constraints on the nHz sGWB provided by the latest PTA measurements can give valuable information about the population of low-mass MBHBs (< 10⁷ M_⊙) that will be detectable by the LISA space-based mission. To this end, we made use of the L-Galaxies semi-analytical model which runs on top of the Millennium suite of simulations and includes detailed physical models to trace galaxy, MBH, and MBHB formation and evolution.

The starting point consisted in creating a population of MBHBs, which produces a nHz sGWB amplitude compatible with the latest PTA results (A = 1.7 − 3.2 × 10⁻¹⁵). To this end, we followed EPTA Collaboration & InPTA Collaboration (2024) which pointed out that the recent sGWB would imply a faster and larger mass growth of MBHs than usually assumed. However, Izquierdo-Villalba et al. (2023) showed that raising the growth efficiency of MBHs to match a louder sGWB tends to over-predict key electromagnetic constraints such as the quasar bolometric luminosity function or the local MBH mass function. To avoid this shortcoming and taking as reference the fiducial MBH growth model of L-Galaxies, we constructed a new framework in which the increase of the galaxy merger efficiency in fuelling gas onto MBHs triggers super-Eddington accretion events. This new model, calibrated by making use of the large volume provided by the Millennium simulation, allowed us to create a population of MBHs which generates a sGWB amplitude of 1.8 × 10⁻¹⁵ and reduces the tension with the electromagnetic constraints on the quasar luminosity function.

With this new model, denoted as Loud, we explored the predictions for the LISA MBHB merger rate and compared them against our former fiducial model (Quiet case), which produces a smaller sGWB amplitude. To reach the range of MBHB masses targeted by LISA, we applied our models on top of the Millennium-II simulation (i.e. the high-resolution version of Millennium) whose merger trees offer the possibility of tracing the cosmological assembly of galaxies and MBHs placed in halos of [10⁷ − 10¹⁴] M_⊙. The main results can be summarised as follows:

• The overall LISA detection rate of MBHBs is not significantly affected by the underlying PTA signal. The Loud model predicts a LISA MBHB detection rate of 12.3 yr⁻¹ whereas the Quiet model forecasts 12.2 yr⁻¹. Therefore, under the assumption of a faster MBH assembly, our results suggest that LISA rates cannot be constrained by using the latest nHz sGWB results.

• The underlying PTA signal causes some differences in the mass distribution of the detected LISA MBHBs. In the Loud model the number of 10^3 − 5 M_⊙ detectable MBHBs decreases by a factor of 1.2 with respect to what is predicted by the Quiet model. Conversely, the number of coalescence of 10^5 − 7 M_⊙ MBHBs is boosted by a factor of 1.5.

• The increase of detectable merging MBHBs with masses > 10⁵ M_⊙ found in the Loud model implies better prospects for multimessenger astronomy. Specifically, the model predicts that MBHBs of 10^5 − 7 M_⊙ potentially detected by LISA (S/N > 10) have 50% (10%) probability of displaying an electromagnetic emission with L_bol > 10⁴³ erg s⁻¹ (L_bol > 10⁴⁴ erg s⁻¹). Furthermore, the LISA detection rate of such type of systems at z < 3 is expected to be 2.4 yr⁻¹.

The results listed above point out that, under the assumption of a faster MBH assembly, the PTA signal cannot constrain the expected LISA merger rate. Indeed, we have shown that a fast MBH growth can only be attained in already-massive systems, where gas is efficiently fuelled onto > 10⁶ M_⊙ MBHs and prompted to lead super-Eddington accretion episodes. On the contrary, MBHs with mass < 10⁶ M_⊙ (i.e. the main contributors to LISA events) are not able to increase their mass as fast as the larger MBHs because the small galaxies where they reside are not capable of sustaining massive and continuous gas inflows towards their centres after galaxy interactions. As a result, the population of low-mass MBHs and MBHBs is just mildly affected by an efficient growth model with episodic super-Eddington accretion events. Consequently, the latter model only leads to small changes in the global LISA merger rate, despite producing a 1.5 louder sGWB at nHz frequencies with respect to our fiducial model. Regardless of these small differences, our results show that an efficient mass-growth model induces an increased number of merging systems with > 10⁶ M_⊙ that can be effectively detected by LISA. Interestingly, MBHs of these masses are also the systems which are more prone to exhibit detectable EM counterparts easily accessible with current and future astronomical facilities. Therefore, our work suggest that, if the astrophysical nature of the PTA signal is confirmed, the possibility of performing multi-messenger analysis with MBHBs could be larger than currently envisioned.

Finally, we stress that in this work we have included a fast assembly of the MBHs to reach the recent sGWB level reported by the PTA collaborations. Nonetheless, different dynamical models for MBHBs could result in a similar enhancement of the nHz signal without invoking any change to the whole MBH population. In an upcoming paper, we will explore this possibility by investigating the conditions leading to a rise in the nHz signal as a consequence of a modified description of the MBHB dynamics.

Acknowledgments

We thank the B-Massive group at Milano-Bicocca University for useful discussions and comments. D.I.V. and A.S. acknowledge the financial support provided under the European Union’s H2020 ERC Consolidator Grant “Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement: 818691). M.C. acknowledges funding from MIUR under the grant PRIN 2017-MB8AEZ and from the INFN TEONGRAV initiative. D.S. acknowledges funding from National Key R&D Program of China (grant no. 2018YFA0404503), the National Science Foundation of China (grant no. 12073014), the science research grants from the China Manned Space project with no. CMS-CSST2021-A05 and Tsinghua University Initiative Scientific Research Program (no. 20223080023). S.B. acknowledges support from the Spanish Ministerio de Ciencia e Innovación through project PID2021-124243NB-C21. M.B. acknowledges support provided by MUR under grant “PNRR – Missione 4 Istruzione e Ricerca – Componente 2 Dalla Ricerca all’Impresa – Investimento 1.2 Finanziamento di progetti presentati da giovani ricercatori ID:SOE_0163” and by University of Milano-Bicocca under grant “2022-NAZ-0482/B”.

References

Appendix A: Effect of the thresholds

In Fig. A.1 we explore how different values of ℛ^th and ${\mathrm{M}}_{\mathrm{inflow}}^{\mathrm{th}}$ affect the number of super-critical accretion events and, thus, the evolution of the luminosity functions. Notice that when varying a threshold we do not impose any limit for the other. In this way, it is possible to marginalize the impact that each parameter. Regarding the effect of ℛ^th, we can see that the smaller the value, the larger the number density of z > 3 AGNs. This is because a small value of ℛ^th permits a large fraction of MBHs that are embedded in relatively gas-poor environments to trigger super-critical accretion. As a result, the population of > 10⁷ M_⊙ MBHs (L_bol > 10⁴⁵ erg/s) is in place earlier in the Universe, having the possibility of triggering more and brighter AGNs than in a case of a model run with large ℛ^th. Interestingly, small values of ℛ^th (e.g., 1) have an opposite effect at low-z, i.e. the normalization of the LFs is smaller than for the cases of large ℛ^th thresholds (e.g., > 10²). This smaller number density is the consequence of the faster assembly of MBHs which consumed most of their gas reservoirs at high-z and thus became inactive (or quiescent) at lower redshifts. Regarding the effect of ${\mathrm{M}}_{\mathrm{inflow}}^{\mathrm{th}}$, we can see similar trends to the ones shown for ℛ^th. Allowing that small inflows fuel super-critical accretion causes a faster assembly of a large fraction of the MBH population, and a rise, at z > 3, in the normalization of the LF at any bolometric luminosity. The drawback of the fast MBH assembly with small ${\mathrm{M}}_{\mathrm{inflow}}^{\mathrm{th}}$ is that at low-z (z < 2) the number of AGNs is diminished, a result of the fact that MBHs consumed more of their reservoirs at high-z.

Fig. A.1. AGN luminosity functions at z = 1.0, 2.0, 4.0 for the model which includes super-Eddington events. The merger trees used correspond to the ones of Millennium. The error bars display the Poissionian error. The results are compared with the observations of Hopkins et al. (2007) (triangles) Aird et al. (2015) (squares) and Shen et al. (2020) (circles). The upper panels correspond to the predictions when varying ℛth and not imposing any ${\mathrm{M}}_{\mathrm{inflow}}^{\mathrm{th}}$ limit. The lower ones represent the LFs of the model when changing ${\mathrm{M}}_{\mathrm{inflow}}^{\mathrm{th}}$ without any ℛth threshold.

Appendix B: The population of massive black holes

In this appendix, we present the population of MBHs generated by L-Galaxies and Millennium merger trees by making use of the Quiet and Loud model. Fig. B.1 depicts the evolution of the black hole mass function. As we can see, regardless of the model, the predictions at z = 0 are consistent with the observational constraints provided by Marconi et al. (2004), Shankar et al. (2004) and Shankar et al. (2013). As described in Izquierdo-Villalba et al. (2020), the mass function displays a fast growth until z ≲ 1 moment at which it slows down and a very small evolution is seen in the massive end (M_BH > 10⁷ M_⊙) and the main evolutionary role is taken by the small MBH population (10^5.5 < M_BH < 10⁷ M_⊙). As expected, the biggest difference between Quiet and Loud is that the latter displays a faster assembly of the massive end of the black hole mass function. Specifically, at z ∼ 5 − 6 the Loud displays a population of MBHs (10^7 − 8 M_⊙) that is absent in the Quiet case.

Fig. B.1. Redshift evolution of the black hole mass function compared to the observational results of Marconi et al. (2004), Shankar et al. (2004, 2009) and Shankar et al. (2013) for the Loud model.

Appendix C: Implications for the MBHs: Rare events in particular hosts

Since Loud model has been calibrated using electromagnetic and GW constraints we can make predictions about the expected comoving number density of MBHs undergoing a Super-Eddington phase at different cosmological times. Fig. C.1 shows the results. As we can see, the vast majority of active MBHs are in a thin disc regime (0.03 < f_Edd < 1), with a peak occurring at 1 < z < 2 coinciding with the peak of star formation and galaxy mergers. The number density of MBHs undergoing a super-Eddington phase is up to one order of magnitude smaller and its shape displays a peak at 3 < z < 4. This maximum is followed by a sharp decrease with a number densities below 10⁻⁵ Mpc⁻³ at z ∼ 2. In the upper panel of Fig. C.1 the population has been divided between massive (> 10⁶ M_⊙) and light (< 10⁶ M_⊙) MBHs. Interestingly, only a small fraction of light MBHs undergo a super-critical accretion. This is presumably caused by the fact that the conditions required to sustain super-Eddington episodes are not easily reachable by the small galaxies where these light MBHs reside. On the contrary, for MBHs of masses > 10⁶ M_⊙ the requirements for super-Eddington episodes are easier to full-field provoking that the relative difference between AGNs triggered by thin discs and super-critical accretion episodes to be smaller than in the case of light MBHs.

Fig. C.1. Comparison between the population of MBHs accreting during the thin disc mode (0.03 < fEdd ≤ 1.0) and super-Eddington (fEdd > 1). The results correspond to the ones of L-Galaxies applied on the Millennium simulation. Upper panel: Number density of AGNs accreting at super-Eddington (red) and at thin disc (black) regimen. While circles represent the whole population of AGNs, squares, and stars correspond to AGNs triggered by MBHs with MBH > 106 M⊙ and MBH < 106 M⊙, respectively. Middle and bottom panels: Median stellar mass (specific star formation rate, sSFR) of the galaxies hosting AGNs in the super-Eddington (red) and thin disc regimen (black). The shaded area corresponds to the percentile 16th–84th.

Regarding the hosts of the MBHs undergoing super-Eddington accretion, Fig. C.1 shows the median stellar mass as a function of redshift. As shown, at high-z (z > 3) the population harboring super-Eddington MBHs does not display differences with respect to the one in a thin disc regime: the hosts display a stellar mass of 10^{8.5 − 9.5} M_⊙. For lower redshifts, the galaxies where super-Eddington MBHs are placed are up to 1 dex more massive than the ones of normal AGNs (10⁹ M_⊙ versus 10^10.5 M_⊙). Besides stellar mass, Fig. C.1 shows the specific star formation rate (sSFR) of galaxies hosting super-Eddington and thin disc accretion. Interestingly, the former displays systematically larger values. This is the result of the fact that large galaxy inflows able to trigger super-critical accretion onto MBHs are principally related to gas-rich major mergers which, in turn, are linked with intense bursts of star formation.

All Tables

All Figures

Fig. 1. Quasar luminosity functions at z = 0.5, 1.0, 2.0, 3.0 for the Quiet (blue) and Boosted model (red) run in the Millennium merger trees (i.e. L-Galaxies + MS). The error bars correspond to the Poissonian error. The results are compared with the observations of Hopkins et al. (2007) (triangles), Aird et al. (2015) (squares) and Shen et al. (2020) (circles). The lower panel represents the ratio between the model and the data points.

In the text

Fig. 2. Minfolw − MRes(t0)/MBH(t0) plane in three different redshifts bins according to L-Galaxies run on the Millennium merger trees. While the y-axis depicts how gas-rich is the environment around the MBH, the x-axis illustrates how powerful are the gas inflows towards the MBH. The blue horizontal and vertical lines correspond to the chosen ℛth, and $M_{\mathrm{inflow}}^{\mathrm{th}}$. Each row corresponds to each of the events that can trigger cold gas inflows toward the MBH: major mergers (top), minor mergers (middle) and disc instabilities (bottom). The plot has been computed by using the free parameters of the Boosted model reported in Table 1. To guide the reader, major (minor) mergers correspond to galaxy interactions which involve galaxies with baryonic masses differing by less (more) than a factor of 2.

In the text

Fig. 3. Evolution of the quasar luminosity function produced by the Loud model (black line), Boosted model (red line) and Quiet model (blue line) when run on the Millennium merger trees. The error bars correspond to the Poissionian error. The results are compared with the observations of Hopkins et al. (2007) (orange triangles), Aird et al. (2015) (purple squares) and Shen et al. (2020) (blue circles). The panels below each luminosity function represent the ratio of the luminosity functions predicted by Loud and Boosted model. Notice that the two lower ones do not display any line because their ratio goes beyond the limit 1.

In the text

Fig. 4. Number density of active and inactive MBHs/MBHBs. Upper panel: number density of MBHs with mass > 108 M⊙ as a function of redshift, z, for Loud (black) and Boosted model (grey). Middle panel: number density of active (fEdd > 0.01) and inactive (fEdd < 0.01) MBHs with mass > 108 M⊙ as a function of redshift, z, for Loud (black) and Boosted model (grey). Lower panel: number density of MBHBs with chirp mass ℳ > 108 M⊙ as a function of redshift, z, for Loud (black) model. For reference, it has been shown the results for the Quiet model (grey). In all panels, the error bars correspond to the Poisson error. We remind the reader that the comoving volume provided by the Millennium simulation corresponds to ∼3.6 × 108 Mpc3.

In the text

Fig. 5. Distribution of MBHB mergers predicted by the model Loud in the plane redshift (z) versus total rest-frame MBHB mass (MBin). Each pixel is encoded by the median signal-to-noise ratio (S/N). The contours represent the number of objects in the z − MBin plane. The small plots around the z − MBin plane correspond to the redshift (right panel) and source-frame mass (top panel) distribution of the differential number of MBHBs. In the right panel, the total distribution has been divided into 5 different mass bins.

In the text

	Fig. 6. Distribution of the z − MBin when each pixel is encoded by the ratio of the number of MBHB mergers in the Loud and Quiet model.
In the text

Fig. 7. Cumulative distribution function (CDF) of bolometric (black) and hard X-ray (2 − 10 KeV, red) luminosity of LISA detectable (S/N > 10) MBHBs with masses MBin > 105 M⊙. Darker lines correspond to the Loud model whereas the light curves depict the predictions of the Quiet one. Horizontal lines highlight the CFR values of 0.5 and 0.2. The vertical pink line highlights the luminosity value of 1043 erg s−1, while black dashed and dotted lines represent the minimum luminosity that a source must have respectively at z = 1 and z = 0.5 to be detected by Athena X-ray observatory assuming a flux limit in the hard X-ray of 2 × 10−15 erg s−1 cm−2 (see Lops et al. 2023).

In the text

Fig. A.1. AGN luminosity functions at z = 1.0, 2.0, 4.0 for the model which includes super-Eddington events. The merger trees used correspond to the ones of Millennium. The error bars display the Poissionian error. The results are compared with the observations of Hopkins et al. (2007) (triangles) Aird et al. (2015) (squares) and Shen et al. (2020) (circles). The upper panels correspond to the predictions when varying ℛth and not imposing any ${\mathrm{M}}_{\mathrm{inflow}}^{\mathrm{th}}$ limit. The lower ones represent the LFs of the model when changing ${\mathrm{M}}_{\mathrm{inflow}}^{\mathrm{th}}$ without any ℛth threshold.

In the text

	Fig. B.1. Redshift evolution of the black hole mass function compared to the observational results of Marconi et al. (2004), Shankar et al. (2004, 2009) and Shankar et al. (2013) for the Loud model.
In the text

Fig. C.1. Comparison between the population of MBHs accreting during the thin disc mode (0.03 < fEdd ≤ 1.0) and super-Eddington (fEdd > 1). The results correspond to the ones of L-Galaxies applied on the Millennium simulation. Upper panel: Number density of AGNs accreting at super-Eddington (red) and at thin disc (black) regimen. While circles represent the whole population of AGNs, squares, and stars correspond to AGNs triggered by MBHs with MBH > 106 M⊙ and MBH < 106 M⊙, respectively. Middle and bottom panels: Median stellar mass (specific star formation rate, sSFR) of the galaxies hosting AGNs in the super-Eddington (red) and thin disc regimen (black). The shaded area corresponds to the percentile 16th–84th.

In the text