© The Authors 2023

1. Introduction

During the last few decades, astronomical observations have proven the existence of massive black holes (MBHs) above 10⁵ M_⊙ lurking at the centre of most of the galaxies (see e.g. Sargent et al. 1978; Tonry 1984; Dressler & Richstone 1988; Kormendy & Richstone 1995; Magorrian et al. 1998; Greene et al. 2020). Furthermore, the masses of these objects show correlations with the properties of their host galaxies, pointing to a co-evolution between MBHs and galaxies (Graham et al. 2001; Häring & Rix 2004; Kormendy & Ho 2013; Savorgnan et al. 2016; Capuzzo-Dolcetta & Tosta e Melo 2017). According to our current hierarchical structure formation paradigm, the assembly of galaxies takes place mainly through a series of mergers and the accumulation of intergalactic gas funnelled into dark matter filaments (White & Rees 1978; White & Frenk 1991; Kauffmann et al. 1999). Therefore, the presence of MBHs in most of the merging galaxies suggests the possible existence of massive black hole binary systems (MBHBs, Begelman et al. 1980). To date, many different processes have been proposed as plausible mechanisms for the formation and evolution of MBHBs. Dynamical friction with background stars, interaction with massive gas clumps, or torques exerted by galactic bars seem to regulate the evolution of MBHs at galactic scales (∼ kpc, Milosavljević & Merritt 2001; Yu 2002; Mayer et al. 2007; Fiacconi et al. 2013; Bortolas et al. 2020, 2022; Li et al. 2022). At smaller distances (< pc), interactions with individual stars, circumbinary gaseous discs, or emission of gravitational waves (GWs) tighten the MBHB and drive the two MBHs to a 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; Bonetti et al. 2020; Franchini et al. 2021, 2022).

According to General Relativity, coalescing MBHBs are among the loudest sources of low-frequency GWs. The future space-based mission LISA (Laser Interferometer Space Antenna, Amaro-Seoane et al. 2017) will tackle the search for GW signals at 0.1 − 100 mHz from coalescing MBHBs of 10⁴ − 10⁷ M_⊙. The population of MBHBs detected by LISA will help in answering several open questions. For instance, the discovery of 10⁴ − 10⁷ M_⊙ binary systems at z < 9 would help us to identify the mechanisms that trigger the growth of low-mass MBHs (Lupi et al. 2016; Pezzulli et al. 2017; Trinca et al. 2022; Spinoso et al. 2023; Sassano et al. 2023). In addition, the detection of MBHs in the first stages of their evolution will shed light on their still unknown origin (Loeb & Rasio 1994; Schneider et al. 2002; Yoshida et al. 2003; Koushiappas et al. 2004; Agarwal et al. 2012; Valiante et al. 2017; Visbal & Haiman 2018; Mayer & Bonoli 2019; Lupi et al. 2021; Volonteri et al. 2021; Spinoso et al. 2023).

Another important science goal of LISA is to use MBHB mergers as ‘bright standard sirens’. The simultaneous observation of LISA sources in the GW and electromagnetic (EM) band will allow the constraint of the luminosity distance–redshift relation and provide an independent measure of the Hubble expansion parameter (Petiteau et al. 2011; Tamanini et al. 2016). In particular, the search for the EM counterpart of the merging MBHB is crucially dependent on LISA’s sky-localisation capabilities. The latest estimates suggest that the sky area constrained by LISA can range from several hundreds of deg² to fractions of deg² depending on the intrinsic properties of the binary, its luminosity distance, and the detection time prior to the final coalescence (Mangiagli et al. 2020; Marsat et al. 2021; Piro et al. 2023). The recent work of Lops et al. (2023) proved that the unequivocal identification of the galaxy housing LISA MBHB will be challenging (see also Kocsis et al. 2006). Specifically, using simulated data, the authors showed that the LISA sky-localisation area is expected to be crowded, with as many as 10⁵ potential host candidates, especially for the fields of z > 1 MBHBs. Lops et al. (2023) highlighted that a pre-selection of X-ray AGNs would reduce this number by one order of magnitude. Still, the large number of galaxies lying within the LISA sky area will hamper the possibility of using LISA MBHBs as standard sirens. Therefore, further theoretical studies are required to provide information as to specific galaxy properties that can be used to unequivocally identify the galaxies housing LISA systems. Several attempts have been made in this direction. For instance, using the cosmological hydrodynamical simulation Illustris, DeGraf et al. (2021) investigated the mergers of MBHs and their connection with the morphologies of the galaxies in which they are found. Assuming an instantaneous MBH coalescence immediately after the galaxy merger, it was shown that the hosts of LISA-like MBHBs will exhibit short-lived merger morphologies (∼500 Myr). Furthermore, these authors showed that with the incorporation of more realistic dynamics of MBHB, the host galaxy has enough time to relax prior to the emission of GW, which means that the connection between coalescing MBHs and post-merger galaxies is blurred even more. Similar results were reported by Volonteri et al. (2020), who analysed the NewHorizon simulation. Moreover, these authors showed that the long time delays between the galaxy and MBHB merger will mean that disturbed features present in the galaxies hosting GW sources will not correlate at all with the merger that led to the MBHB formation.

Under these premises, in this work, we explore the properties of the galaxies hosting z < 3 LISA MBHBs and use a statistical study to investigate whether or not signatures of galaxy mergers can be used as tracers for seeking the hosts of merging MBHBs. To this end, we use the L-Galaxies semi-analytical model (SAM, Henriques et al. 2015; Izquierdo-Villalba et al. 2020, 2022a; Spinoso et al. 2023) applied to high-resolution dark matter merger trees of the Millennium-II simulation (Boylan-Kolchin et al. 2009) to produce a synthetic lightcone of ∼1000 deg² reaching z ∼ 3.5. Thanks to the physics included in L-Galaxies, the resulting lightcone is detailed enough to account for the cosmological evolution of galaxies and of their single and binary MBHs. We highlight that throughout the whole paper, the LISA MBHB population is chosen according to the total mass of the binary, its mass ratio, and its time to coalescence. Specifically, we impose the total mass 10^4 − 7 M_⊙, a mass ratio of ≥0.1, and a merging time of < 1 Myr to ensure that the selected MBHBs eventually emit GW signals inside the LISA frequency band.

The paper is organised as follows: In Sect. 2 we describe the main characteristics of the Millennium-II dark matter simulation and summarise the physics implemented in L-Galaxies to tackle the assembly and evolution of galaxies, MBHs, and MBHBs. In Sect. 3 we present the abundance and properties of MBHBs potentially detectable by the LISA interferometer. In Sect. 4 we summarise the properties of the galaxies hosting LISA MBHBs and compare them with those of galaxies housing single MBHs with the same mass as the binary systems. In Sect. 5 we discuss the possibility of using merger signatures as a potential distinctive feature of LISA hosts; that is, the galaxies where detectable LISA MBHBs are found. In Sect. 6 we list several caveats that should be taken into account when interpreting our results. Finally, in Sect. 7 we summarise our main findings. A Λ cold dark matter (ΛCDM) cosmology with parameters Ω_m = 0.315, Ω_Λ = 0.685, Ω_b = 0.045, σ₈ = 0.9, and H₀ = 67.3 km s⁻¹ Mpc⁻¹ is adopted throughout the paper (Planck Collaboration XVI 2014).

2. A lightcone for the study of the LISA massive black hole binaries

In this section, we describe the dark matter (DM) simulation and galaxy formation model used to generate a lightcone specifically designed to study potential LISA MBHB sources. We use the so-called L-Galaxies semi-analytical model (SAM), a state-of-the-art model set to reproduce many different observational constraints such as the stellar mass function, the cosmic star formation rate density evolution, galaxy colors, and the fraction of passive galaxies (we refer to Guo et al. 2011; Henriques et al. 2015 for further details). Among all the versions of the model, we use the one presented in Henriques et al. (2015) with the modifications of Izquierdo-Villalba et al. (2019a, 2020, 2022a). These changes were included to improve the predictions for galaxy morphology, extend the physics of MBHs, and introduce the formation and evolution of MBHBs. In the following, we summarise the main features of the model, and we refer the reader to the papers cited above for a detailed description of the baryonic physics included.

2.1. The underlying dark matter population: Millennium-II

L-Galaxies is a semi-analytical model, which self-consistently couples different astrophysical processes with the DM merger trees of N-body simulations. In particular, our SAM can be run on top of the merger trees of the Millennium suite of simulations: Millennium-I (MS, Springel 2005), Millennium-II (MSII, Boylan-Kolchin et al. 2006), and Millennium-XXL (MXXL, Angulo et al. 2012). The different box sizes and DM mass resolution of the Millennium suite offer the possibility to explore different baryonic processes over a wide range of scales and environments.

Among all these Millennium simulations, we use the MSII, whose mass resolution allows the study of MBHs and MBHBs hosted in galaxies with stellar mass as low as ∼ 10⁷ M_⊙. Briefly, the MSII follows the cosmological evolution of 2160³ DM particles with mass 6.885 × 10⁶ M_⊙ h⁻¹ within a periodic comoving box of 100 M_⊙ h⁻¹ on a side. The simulation was stored at 68 different epochs or ‘snapshots’, to which the SUBFIND algorithm was applied to detect all the DM halos whose minimum halo mass corresponds to seven times the particle mass (∼10⁷ M_⊙ h⁻¹). Subsequently, the L-HALOTREE code arranged these structures according to their evolutionary path in the so-called merger trees (Springel et al. 2001; Springel 2005). The time resolution offered by the finite number of the MSII outputs leads to inconveniences in accurately tracing the baryonic physics involved in galaxy evolution. To overcome this, the SAM does an internal time interpolation between two consecutive snapshots with ∼5 − 20 Myr of time resolution depending on redshift. Finally, L-Galaxies rescales the original cosmology of MSII (WMAP1 and 2dFGRS ‘concordance’ ΛCDM framework; Spergel et al. 2003) to the cosmology of Planck first-year data release (Planck Collaboration XVI 2014) using the Angulo & White (2010) methodology. This rescaling modifies the MSII box size and particle mass by factors of 0.96 and 1.12, respectively. Taking into account this, the merger trees of MSII enables us to trace the cosmological assembly of galaxies placed in halos of 5.7 × 10⁷ − 3 × 10¹⁴ M_⊙.

2.2. The assembly of the galaxy population: L-Galaxies formation and evolution model

Following the framework of White & Frenk (1991), L-Galaxies assumes that the birth of a galaxy takes place at the centre of every newly formed DM halo. As soon as a DM halo collapses, a fraction of baryonic matter (proportional to the baryon fraction) is trapped and collapses with it. During this process, the material is shock heated and forms a diffuse, spherical, and quasi-static hot gas atmosphere with an extension equal to the halo virial radius (R_vir). Part of this hot gas is allowed to cool down and migrate towards the DM halo centre (White & Rees 1978). The rate at which this process takes place is determined by cooling functions (Sutherland & Dopita 1993) and the amount of hot gas enclosed within the halo cooling radius (r_cool), which is defined as the radius at which the cooling time matches the halo dynamical time. At high-z and in low-mass DM halos, the hot gas can cool rapidly (rapid infall, r_cool > R_vir), causing the migration of the whole mass towards the DM halo centre at approximately the free-fall rate. On the other hand, a slow ‘cooling flow regime’ (r_cool < R_vir) takes place at low-z and in massive halos. In these cases, only a fraction of the hot gas is allowed to condense through cooling flows. After any of these condensation processes, the cold gas settles in a disc with a specific angular momentum inherited from the host DM halo. This newly formed disc is assumed to be distributed with an exponential profile, whose extension is determined according to the evolution of the gas angular momentum (see Guo et al. 2011).

Based on the observational results of Kennicutt (1998), our SAM assumes that star formation processes take place as soon as the surface density of the cold gas exceeds a critical value. When this occurs, the galaxy begins (or continues) the assembly of its stellar disc on a timescale given by the cold gas disc dynamical time. As a consequence of star formation processes, massive and short-lived stars explode as supernovae (SNe), injecting energy and metals into the cold gas disc (SN feedback). This injection causes the reheating of a fraction of cold gas and may additionally expel a fraction of the hot gas beyond the halo virial radius. At later times, this ejected gas can be reincorporated, initiating new star formation events. In addition to SNe feedback, L-Galaxies introduces the so-called radio-mode feedback as an additional process to regulate the assembly of the stellar component in massive galaxies. This mechanism is activated by the gas accretion onto the MBH from the hot gas atmosphere around the galaxy. The result of this accretion is the release of kinetic energy, whose injection into the surrounding medium can reduce or even suppress the cooling of gas. L-Galaxies computes the extent of the stellar disc by following the evolution of its specific angular momentum (modified by star formation events and galaxy mergers) and assuming an exponential profile (see Guo et al. 2011). Finally, our SAM models large-scale effects (also known as environmental processes), such as ram pressure stripping or galaxy tidal disruption. These processes occur when the halo hosting the galaxy falls into a larger system. As a consequence, the galaxy can lose its entire hot gas atmosphere and can be deprived of its cold gas and stellar component through tidal forces.

2.3. The assembly of bulges

The continuous assembly of the stellar disc causes some galaxies to undergo disc instability (DI) processes. These events refer to the mechanism by which the stellar disc becomes massive enough to suffer non-axisymmetric instabilities. The eventual result is the formation of a central ellipsoidal component via buckling of the nuclear stellar orbits (Mo et al. 1998). The occurrence of DIs is modelled according to the Efstathiou et al. (1982) analytic stability test (tested against cosmological simulations; see Izquierdo-Villalba et al. 2022b but see also Romeo et al. 2023). When the criterion is satisfied, our SAM triggers the formation (or growth) of a bulge by transferring from the disc the minimum stellar mass needed to make it marginally stable again (i.e. $\mathrm{\Delta }{M}_{\mathrm{stars}}^{\mathrm{DI}}$). The effective radius of the bulge after any DI events is determined by assuming that the transferred stellar mass comes from the innermost part of the disc (see Guo et al. 2011, for further information).

The hierarchical growth of the DM halos also shapes galaxy properties. The interaction between galaxies is ruled by the merger of the parent DM halos. As soon as two DM halos merge, their galaxies merge as well on a timescale given by the dynamical friction presented in Binney & Tremaine (1987). According to the baryonic (stars plus gas) merger ratio of the two interacting galaxies (m_R ≤ 1), L-Galaxies differentiates between major (m_R > 0.2) and minor interactions (m_R < 0.2). Major mergers completely destroy the discs of the two galaxies, giving rise to a pure spheroidal remnant that undergoes a ‘collisional starburst’. Conversely, during minor interactions, the disc of the larger galaxy survives and experiences a burst of star formation, while its bulge integrates the entire stellar mass of the satellite. After any type of merger, L-Galaxies determines the effective radius of the remnant bulge by assuming the conservation of the binding and orbital energy. In addition to these two merger treatments, the model used in this work includes the prescription of ‘smooth accretion’ presented by Izquierdo-Villalba et al. (2019a) in order to deal with the physics of extremely minor mergers. Indeed, these latter authors showed that the inclusion of these processes is important in order to recover the observed morphology of dwarf galaxies (M_stellar ≤ 10⁹ M_⊙) in L-Galaxies when this model is run on top of the MSII simulation. Specifically, smooth accretions take place in satellite galaxies with low binding energy. As a result, the stellar component of the satellite (i.e the bulge plus disc) gets diluted inside the disc of the central galaxy before being able to reach the nucleus. Consequently, the central galaxy loses the possibility of forming (growing) the bulge component.

Based on the ratio between the bulge and total stellar component (known as the bulge-to-total ratio, B/T), L-Galaxies divides the galaxy population into several morphological types:

• (a)

  Ellipticals: Galaxies with B/T > 0.7.

• (b)

  Discs or spirals: Galaxies with B/T < 0.7. The bulge component of these systems is split into three types:

  • Pseudobulges: Galaxies with 0.01 < B/T < 0.7 where more than two-thirds of the bulge mass has been accumulated through disc instabilities.

  • Classical bulges: Galaxies with 0.01 < B/T < 0.7 where less than two-thirds of the bulge mass has been accumulated through disc instabilities. Thus, (minor/major) mergers are the main mechanisms that brought mass into the bulge.

  • Extreme late-type: Bulgeless galaxies or galaxies with B/T < 0.01 (regardless of the process that gave rise to the bulge component).

2.4. The population of massive black holes: seeds, growth, and spin

Thanks to the modifications of Izquierdo-Villalba et al. (2020, 2022a) and Spinoso et al. (2023), L-Galaxies is able to track the formation, growth, and spin evolution of MBHs self-consistently. In the following subsections, we summarise the main physics included in our SAM to model the evolution of MBHs.

2.4.1. The seeding of massive black holes cosmic dawn

The formation of MBH seeds in L-Galaxies has been extensively explored by Spinoso et al. (2023). In particular, the genesis of light seeds (PopIII remnants, Bromm & Larson 2004) was accounted for using a subgrid approach, while the formation of massive seeds (i.e. intermediate-mass and heavy BH-seeds) was addressed by taking into account the spatial variations of the IGM metallicity and the UV-background produced by star formation events (see Inayoshi et al. 2020, for a recent review). With these models applied on L-Galaxies and MSII, the authors showed that the formation of BHs is strongly inhibited at z ≲ 6 − 7 due to the progress of IGM chemical enrichment. In addition, the occupation fraction of newly formed BH seeds showed a dependence on the halo mass and redshift, which was shown to be almost null in halos of < 10⁸ M_⊙ at z < 9. Following these results, in the present work we include a simple empirical BH seeding model where MBHs only form within newly resolved galaxies at z ≥ 7 with a probability, 𝒫, given by:

$$\begin{array}{c}\hfill \mathcal{P}=\mathcal{A}{(1+z)}^{\gamma }\left(\frac{M_{\mathrm{halo}}}{M_{\mathrm{halo}}^{\mathrm{th}}}\right),\end{array}$$(1)

where 𝒜 = 0.015, γ = 7/2, and $M_{\mathrm{halo}}^{\mathrm{th}}=7\times {10}^{10}{M}_{\odot }$ (see also Eq. (9) of Spinoso et al. 2023). To guide the reader, in Fig. 1 we show how 𝒫 varies with redshift and halo mass, M_halo. As shown, at very high-z, the seeding process occurs mainly in low-mass halos. As the redshift decreases the seeding events shift towards higher mass halos. This evolution is assumed to rise as a combination of an early formation of MBHs after the explosion of PopIII stars and a later creation of MBHs via stellar runaway mergers and a direct collapse of pristine gas clouds. Based on Eq. (1), every time that a z ≥ 7 galaxy is formed, we compute the value of 𝒫 and draw a random value ℛ ∈ [0 − 1]. If ℛ > 𝒫, a MBH is placed at the centre of the galaxy whose mass and spin is randomly extracted from 10² − 10⁴ M_⊙ and 0 − 0.998, respectively (see Fig. 3 of Spinoso et al. 2023, for further information about the seed mass choice). To show that this chosen seeding procedure retrieves a reasonable MBH population, in Fig. 2 we present the predicted occupation fraction at z = 0. As shown, the occupation fraction decreases towards low stellar masses, which is consistent with the constraints of Miller et al. (2015). While at M_stellar > 10¹⁰ M_⊙, almost all galaxies host a nuclear MBH, and at M_stellar < 10⁹ M_⊙ less than 80% have one. We stress that we plan to include the full physical seeding model developed by Spinoso et al. (2023) in future work.

Fig. 1. Probability 𝒫 that a given DM halo of mass Mhalo would be seeded by a MBH. Each color corresponds to a different redshift. Vertical and horizontal grey dashed lines highlight the halo mass resolution of MSII and the value 𝒫 = 0.5, respectively.

Fig. 2. Massive black hole occupation fraction at z = 0 as a function of stellar mass (Mstellar). The shaded grey region corresponds to the constraints presented in the observational work of Miller et al. (2015).

2.4.2. The growth of massive black holes

In the semi-analytical model, MBHs can grow via three different channels: cold gas accretion, hot gas accretion, and mergers with other MBHs. In the following lines, we summarise the main assumptions related to gas consumption processes:

• Hot gas accretion: This channel of growth is triggered by continuous gas accretion from the hot gas atmosphere that surrounds the galaxy hosting the MBH (Croton et al. 2006). The rate of accretion is usually orders-of-magnitude below the Eddington limit and is linked with the so-called radio mode feedback, which injects energy into the hot atmosphere, halting the cooling gas inflows that supply gas to the galaxy. In the model, the accretion due to hot gas is determined as (Henriques et al. 2015):

  $$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{BH}}=k_{\mathrm{AGN}}\left(\frac{M_{\mathrm{hot}}}{{10}^{11}{M}_{\odot }}\right)\left(\frac{M_{\mathrm{BH}}}{{10}^8{M}_{\odot }}\right),\end{array}$$(2)

  where M_hot is the total mass of hot gas surrounding the galaxy and k_AGN is a free parameter set to 9 × 10⁻⁵ M_⊙ yr⁻¹ to reproduce the turnover at the massive end of the galaxy stellar mass function.

• Cold gas accretion: This is the main channel driving the black hole growth and is triggered by both galaxy mergers/smooth accretion and disc instability events. In particular, after a galaxy merger or smooth accretion, the nuclear MBH can accumulate a fraction of the galaxy cold gas given by

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

  where z_merger is the redshift of the galaxy merger, M_gas the cold gas mass of the galaxy, V₂₀₀ the virial velocity of the host DM halo, and V_BH, $f_{\mathrm{BH}}^{\mathrm{merger}}$ are two adjustable parameters set to 280 km s⁻¹ and 0.014, respectively. On the other hand, after a disc instability, the black hole accretes an amount of cold gas proportional to the mass of stars that triggered the stellar instability, $\mathrm{\Delta }{M}_{\mathrm{stars}}^{\mathrm{DI}}$:

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

  where z_DI is the redshift at which the disc instability occurs, and $f_{\mathrm{BH}}^{\mathrm{DI}}$ is a free parameter that takes into account the gas-accretion efficiency, and is set to 0.0014. All these adjustable parameters have been tuned to give the best agreement between the observations and model predictions for the z = 0 black hole-bulge correlation and the BH mass function for M_BH > 10⁶ M_⊙. After a galaxy merger or a disc instability, the cold gas available for accretion (see Eqs. (3) and (4)) is assumed to settle in a reservoir around the black hole (with total mass M_Res), which is progressively consumed according to a two-phase model extensively used in Izquierdo-Villalba et al. (2020, 2022a). The first phase corresponds to an Eddington-limited growth, which lasts until the MBH consumes a faction ℱ of the available gas reservoir. The free parameter ℱ is set to 0.7 in order to match the faint end of the low-z AGN LFs (see also Marulli et al. 2008; Bonoli et al. 2009). Once this phase ends, the BH enters a self-regulated or quiescent growth regime characterised by progressively lower accretion rates. To show that the model of MBHs explained in this section gives rise to a population in good agreement with the observations, in Fig. 3 we present the black hole mass function (BHMF) and the bulge–MBH mass correlation in the local Universe. As we can see, the scaling relation raised in our SAM is consistent with the results reported by Kormendy & Ho (2013). Also, the BHMF is in good agreement with the observations at M_BH > 10⁷ M_⊙. For lower masses, L-Galaxies applied in the MSII merger trees finds a steeper increase than the one seen by observations. Interestingly, this rise is also seen in cosmological hydrodynamical simulations such as EAGLE (Rosas-Guevara et al. 2016) or TNG (Habouzit et al. 2021). The large number density of low-mass black holes reported in Fig. 3 has an impact on the evolution of active MBHs. Figure 4 shows the evolution of the bolometric luminosity function (LF) predicted by L-Galaxies. As shown, our SAM applied on the Millennium-II merger trees is able to reproduce the observed trends of bright objects (L_bol ≥ 10⁴⁶ erg s⁻¹). However, it overpredicts the faint end of the LFs (see similar trends in Sijacki et al. 2015; Griffin et al. 2019; Weinberger et al. 2018; Marshall et al. 2020; Trinca et al. 2022). Interestingly, this excess is not seen in runs with the Millennium simulation, pointing out that the evolution of the faint end of the LFm, and therefore the low-mass MBH population, is strongly affected by the resolution of the underlying DM simulation. The nature of the inconsistency between simulated and observed LFs is still an open issue. On one hand, it is challenging to cover wide sky areas with large depths, and therefore a robust sampling of high-z faint AGNs is still missing (Siana et al. 2008; Masters et al. 2012; McGreer et al. 2013; Niida et al. 2016; Akiyama et al. 2018). On the other hand, some studies have explored different ways to suppress the large excess of faint AGNs seen in most of the SAMs and hydro-dynamical simulations (see e.g. Hirschmann et al. 2014; Griffin et al. 2019; Habouzit et al. 2022); for instance, varying the efficiency of MBH seed formation or the ability of newly born MBH to accrete matter have been proposed as plausible means to this end (see e.g. Degraf et al. 2010; Fanidakis et al. 2012; DeGraf & Sijacki 2020; Spinoso et al. 2023; Trinca et al. 2022). Nevertheless, no clear answer has yet been found, and further investigation is needed.

Fig. 3. General properties of MBHs in L-Galaxies SAM applied in the Millennium-II merger trees. Left panel: 2D histogram of the bulge-black hole correlation at z = 0. The predictions are compared with the observations of Kormendy & Ho (2013; K&H2013) where blue, orange, and red points correspond to pseudobulges, classical bulges, and ellipticals, respectively. Right panel: black hole mass function at z = 0 when the model is run on top of the Millennium-II merger trees (black lines). Error bars correspond to the Poissonian error. For comparison, we have added the observational constraints of Marconi et al. (2004) and Shankar et al. (2004, 2009, 2013).

Fig. 4. Quasar bolometric luminosity functions (Lbol) at z ∼ 0.5, 1.0, 2.0, and 3.0 predicted by L-Galaxies when MS (red) and MSII (black) DM merger trees are used. For comparison, the grey dashed lines represent the luminosity functions of L-Galaxies applied in addition to MSII when a cut of ≥107 M⊙ in black hole mass is performed. For comparison, these predictions are compared with the data presented in Shen et al. (2020; blue circles). In all the plots, error bars correspond to the Poissonian error (large error bars correspond to the cases where only one object is found).

2.5. The population of massive black hole binaries

Besides including a comprehensive model for MBH growth, L-Galaxies deals with the dynamical evolution and growth of MBH binaries (MBHBs; see Izquierdo-Villalba et al. 2022a). In the following sections, we summarise key aspects of the model included in L-Galaxies.

2.5.1. The dynamical evolution

The evolution of MBHBs inside L-Galaxies is divided into three different stages (Begelman et al. 1980): pairing, hardening, and the gravitational wave phase.

• Pairing phase: After a galaxy merger, the MBH hosted by the satellite galaxy experiences dynamical friction that leads to its pairing with the nuclear MBH of the primary galaxy. This process causes the satellite MBH (typically deposited at a few kiloparsec) to sink towards the galactic centre of its new host. The time needed for the satellite MBH to reach the nuclear part of the galaxy, $t_{\mathrm{dyn}}^{\mathrm{BH}}$, is determined using (Binney & Tremaine 2008):

  $$\begin{array}{c}\hfill t_{\mathrm{dyn}}^{\mathrm{BH}}=19f_{\mathrm{s}}f(\epsilon ){\left(\frac{r_0}{5\mathrm{kpc}}\right)}^2\left(\frac{\sigma }{200\mathrm{km}{\mathrm{s}}^{-1}}\right)\left(\frac{{10}^8{M}_{\odot }}{M_{\mathrm{BH}}}\right)\frac{1}{\mathrm{\Lambda }}[\mathrm{Gyr}],\end{array}$$(5)

  where f(ε) = ε^0.78 is a function that depends on the orbital circularity of the MBH (ε, Lacey & Cole 1993), r₀ is the initial position of the black hole deposited by the satellite galaxy after the merger, σ is the velocity dispersion of the remnant galaxy, M_BH is the mass of the satellite black hole, and Λ = ln(1 + M_stellar/M_BH) is the Coulomb logarithm (Mo et al. 2010). The computation of all these quantities can be found in Izquierdo-Villalba et al. (2022a). The variable f_s takes into account the stochastic inspiralling of MBHs seen in simulations of clumpy (gas-rich) and barred galaxies. For instance, the simulations of gas-rich galaxies of Tamburello et al. (2017) showed that the interaction between MBHs and massive clumps typically lags the pairing phase of MBHs. At the other extreme, we find the results of Lupi et al. (2015), who by simulating the late stages of a gas-rich galaxy merger found that the gravitational torques after the interaction can be very efficient in forming massive gas clumps that substantially perturb the orbits of the infalling MBHs. These perturbations result in impulsive kicks that lead to the formation of a gravitationally bound MBHB in ∼10 Myr from the start of the merger. Concerning galactic structures, the works of Bortolas et al. (2020, 2022) showed that bar structures induce an erratic motion in pairing MBHs, causing either a delay or a boost to the inspiral. Taking these studies into account, we assume that f_s is set to 1 when the galaxy is either gas poor (f_gas < 0.5) or does not display a pseudobulge structure (i.e. bar-related morphology). For the other cases, f_s is a random value extracted from a log-normal distribution whose free parameters are a median of 0.2 and a variance of 0.6. The choice of these values is motivated by the shape of the resulting log-normal distribution, which peaks at ∼1 and features a positive skewness (i.e. a long tail towards values > 1). We stress that during this pairing phase, the two MBHs do not form a bound system. Instead, they can be considered as dual MBHs or dual AGNs in cases where both of them undergo an active phase.

• Hardening and gravitational phase: Once the dynamical friction phase ends, the satellite MBH reaches the galactic nucleus of the new galaxy and binds with the central MBH (∼pc separation) forming a MBHB. Hereafter, we refer to the most massive black hole in the system as the primary black hole (with mass M_BH, 1), whereas the less massive one is tagged as the secondary black hole (with mass M_BH, 2). The total mass of the binary and its mass ratio are denoted M_Bin = M_BH, 1 + M_BH, 2 and q = M_BH, 2/M_BH, 1, respectively. Regarding the initial properties of the binary orbit, the code assumes that the eccentricity of the binary, e₀, starts with a random value in the range [0 − 1] while the initial separation, a₀, is set to the scale in which M_Bulge(< a₀) = 2 M_BH, 2, where M_Bulge(< a₀) corresponds to the mass in stars of the hosting bulge within a₀. To determine a₀, L-Galaxies assumes that the bulge mass density profile follows a Sérsic model (Sersic 1968) with an index extracted randomly according to the observed distribution of z = 0 pseudobulges, classical bulges, and ellipticals (Gadotti 2009). We stress that no assumptions are needed to determine the normalisation and scale radius of the Sérsic profile because L-Galaxies computes the redshift evolution of the bulge mass and effective radius self-consistently (see Guo et al. 2011; Izquierdo-Villalba et al. 2019a). Once the two MBHs bind at the galactic centre, the separation (a_BH) and eccentricity (e_BH) of the binary system are evolved depending on the environment in which the binary is embedded. If the gas reservoir around the binary (M_Res) is larger than its total (M_Bin), the evolution of the system is driven by the interaction with a circumbinary gaseous disc and then GW emission (Dotti et al. 2015). Otherwise, the system evolves thanks to the interaction with single stars embedded in a Sérsic profile and the emission of GWs (Quinlan & Hernquist 1997; Sesana & Khan 2015). We refer the reader to Izquierdo-Villalba et al. (2022a) for a full description of the equations used to evolve a_BH and e_BH.

2.5.2. Triple interactions

In some cases, the lifetime of a MBHB can be long enough that a third MBH can reach the galaxy nucleus and interact with the binary system. In this scenario, multiple outcomes are allowed. To deal with these triple MBH interactions, L-Galaxies uses the tabulated values of Bonetti et al. (2018). Based on the mass of the intruder MBH and the mass ratio of the MBHBs, the model determines whether the triple interaction leads to the prompt merger of the MBHB or causes the ejection of the lightest MBH from the system. In case the latter scenario takes place, the separation of the leftover MBHB is computed following Volonteri et al. (2003) and the resulting e_BH is selected as a random value in the range [0−1].

2.5.3. The growth of massive black hole binaries

We model the growth of MBHBs in a different way from that of single MBHs. In particular, our SAM assumes that the accretion rates of the two MBHs of the binary are correlated, as proposed by Duffell et al. (2020):

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

where ${\dot{M}}_{{\mathrm{BH}}_1}$ and ${\dot{M}}_{{\mathrm{BH}}_2}$ are respectively the accretion rate of the primary and secondary MBHs. For simplicity, the latter is set to the Eddington limit, as in Izquierdo-Villalba et al. (2022a).

2.5.4. The growth of massive black holes in the pairing phase

On top of the growth of MBHBs, L-Galaxies deals with the gas accretion of MBHs in the dynamical friction phase. For these objects, the gas consumption is modelled in the same way as for nuclear black holes (Sect. 2.4.2). The growth lasts until the MBH consumes the total gas reservoir stored prior to the merger. This reservoir is set as the sum of all the gas that the MBH accumulated before the galaxy merger (i.e. as a consequence of past disc instabilities or mergers) and an extra amount computed at the time of the galaxy merger as Eq. (3), where the cold gas in the equation is assumed to be that of the satellite galaxy. This extra accumulation of gas is motivated by the hydrodynamical simulations of merging galaxies with central MBHs by Capelo et al. (2015). The authors showed that during the merging process, the secondary galaxy suffers important perturbations during the pericentre passages around the central one. In such cases, the black hole of the secondary galaxy experiences accretion enhancements mainly correlated with the galaxy mass ratio.

2.6. Lightcone construction

In this work, we explore the properties of LISA MBHBs and their host galaxies by making use of a simulated lightcone, that is, a mock universe in which only galaxies whose light has just enough time to reach the observer are included. The main limitation to creating a lightcone using L-Galaxies and MSII is the small box-side length of the latter (L ∼ 100 Mpc h⁻¹), which is insufficient to represent the Universe beyond redshift 0.025. To overcome this, we employ the methodology presented in Izquierdo-Villalba et al. (2019b). Briefly, to reach a desired redshift depth, the procedure exploits the periodic boundaries of the MSII and replicates its simulated box (i.e. fundamental box) 𝒩 times in each Cartesian coordinate. Once the replication is made, the methodology establishes the location and line-of-sight (LOS) of the observer. Specifically, the observer is placed at the origin of the first replica, while the LOS is chosen according to Kitzbichler & White (2007) to minimise the structure replication inside the lightcone. As shown by Kitzbichler & White (2007), selecting a LOS given by $\widehat{u}=(n,m,nm)/|(n,m,nm)|$ (where n and m are two different integers with no common factor) implies that the observer will pass through the first periodic image at a distance given by mnL. Furthermore, no point of the fundamental box is imaged more than once within mnL when a rectangular footprint of size 1/m²n × 1/n²m (radians) is selected. Once the position and LOS of the observer are set, the methodology presented in Izquierdo-Villalba et al. (2019b) places galaxies inside the lightcone by determining the moment at which they (and their corresponding single and binary MBHs) cross the past lightcone of the observer. To this end, the galaxy merger trees provided by L-Galaxies were used given that they accurately follow in time the cosmological evolution of individual galaxies between the DM snapshots with a fine time step resolution (see Sect. 2.1).

As we are interested in z < 3 galaxies, we set 𝒩 = 48, n = 5, and m = 9. This selection implies a LOS of (RA, Dec) = (77.1, 60.95) deg. We highlight that with this setup, no structure repetition would be allowed up to z ∼ 3 for a FOV of 0.14 × 0.25 deg². As in this work a large area is required in order to have sufficient statistics, we allow some structure repetition and we set the lightcone footprint as a rectangular shape with an extension of (δRA, δDec) = (45.6, 22.5) deg (corresponding to 1027 deg²).

3. From gravitationally bound systems to potential LISA sources

In this section, we explore the abundance of MBHBs at different masses and distances. Specifically, we define LISA MBHBs (or LISA systems) as those binaries whose total mass is 10⁴ ≤ M_Bin ≤ 10⁷ M_⊙, q ≥ 0.1¹ and whose time to coalescence is < 1 Myr (i.e. GW-dominated phase). These cuts ensure that the selected population emits GWs at 0.1 − 100 mHz.

In Fig. 5 we present the number of MBHBs per deg² as a function of redshift. The population of MBHBs is divided into three different mass bins. The lightest systems, binaries with M_Bin = 10⁴ − 10⁵ M_⊙, are less numerous towards low z. For instance, at z ∼ 3, the number of hard binaries can reach up to 10 deg⁻² while at z ∼ 0.5 it drops down to 1 deg⁻². A similar trend is shown by single MBHs with the same mass but their abundances can be up to about five times higher, regardless of redshift. When the MBHB sample is divided into bins according to the semi-major axis, the large majority of the systems are at a_BHB < 1 pc with very few of them (a number 100 times smaller) at < 10⁻⁴ pc. Specifically, the population of LISA MBHBs coincides with binaries at < 10⁻⁵ pc whose number density does not overpass ∼10⁻² deg⁻² and represent only 1% of the whole population of 10⁴ − 10⁵ M_⊙ hard binaries.

Fig. 5. Number of hard MBHBs (grey solid) and single MBHs (grey dashed) per deg2 (dN(z)/dΩ, black) as a function of redshift. Different colors represent the same but when the population of MBHBs is divided by semi-major axis (aBH). Each panel corresponds to binaries with different total masses: 104 − 105 M⊙ (top), 105 − 106 M⊙ (middle), and 106 − 107 M⊙ (bottom). Overall, the figure shows that LISA MBHBs represent 0.1%−1% of the MBHs between 104 and 107 M⊙.

MBHBs with a total mass of between 10⁵ and 10⁶ M_⊙ can be up to three times more abundant than MHMBs in the previous mass range, with values reaching ∼10 − 30 deg⁻² regardless of redshift. Despite these large numbers, single MBHs with the same mass are a factor 4 more numerous, with approximately 100 objects per deg². Concerning the MBHB separation, the bulk of the population displays a semi-major axis of 0.1 < a_BHB < 1 pc. LISA systems are a much tighter sample with separations of < 10⁻⁴ pc and an abundance that decreases towards low z. For instance, at z ∼ 3 the number of objects per deg² reaches 0.5 while at z ∼ 0.5 it drops down to 0.01. Finally, MBHBs of 10⁶ − 10⁷ M_⊙ display a different trend with respect to the two previous mass bins. The number of objects in the sky rises from z ∼ 3 (10 deg⁻²) down to z ∼ 0.5 (100 deg⁻²), redshift at which the abundance of MBHBs reaches its maximum. Interestingly, single and binary MBHBs show the same abundance at z < 1, pointing out that half of the low-z MBH population members of 10⁶ − 10⁷ M_⊙ reside in binary systems. Finally, at these masses, LISA MBHBs have separations of a_BHB < 10⁻³ pc with an abundance of ∼0.01 deg⁻², independently of redshift.

Overall, the numbers shown in this section point out that ∼20% of low-mass MBHs are expected to be in relatively wide binaries (separated by one or more parsecs). Furthermore, this fraction can rise up to ∼50% for MBHs of 10⁶ − 10⁷ M_⊙ at low z. These predictions are larger than the ones presented in other works. For instance, the study of Volonteri et al. (2003) showed that about 5%−10% of z ∼ 0 halos host an MBHB². Among these, approximately 60% have a separation of > 0.1 kpc and only ∼10% feature an advanced hardening stage (a_BH < 10 pc). Regardless of these differences (most likely caused by the different approaches and assumptions), previous and current studies highlight that MBHBs at parsec scales located in low-mass galaxies make a non-negligible contribution to the MBH population. Despite this, the faint nature and the small separation of these MBHBs will challenge their discovery. Detecting and characterising the presence of Doppler shifting in broad AGN emission lines can be a good avenue to explore in order to detect (sub)parsec MBHBs with current spectroscopy facilities such as SDSS (see e.g. Bogdanović et al. 2009; Tsalmantza et al. 2011; Montuori et al. 2011; Eracleous et al. 2012; Shen et al. 2013). On the other hand, upcoming deep surveys such as LSST (Ivezić et al. 2019) or Athena (Nandra et al. 2013) will help in building a complete census of active low-mass MBHs at cosmological distances. Furthermore, the possibility of getting periodic light curves from the data of these observatories will open a new path towards identifying and characterising potential low-mass MBHB candidates (see e.g. Valtonen et al. 2008; Graham et al. 2015; Charisi et al. 2016; Liu et al. 2016, 2019; Liao et al. 2021; Witt et al. 2022).

4. The hosts of LISA massive black hole binaries

The detection of GWs coming from 10⁴ − 10⁷ M_⊙ merging MBHBs, together with the identification of their host galaxies, will create a new opportunity to study the population of MBHs and constrain our standard cosmological model. In this section, we aim to guide the future search for z < 3 LISA hosts by determining their masses and properties. To this end, in Fig. 6 we present the median stellar mass of the galaxies where LISA MBHBs are placed. As shown, these GW sources inhabit low-mass galaxies (i.e. dwarf range) of M_stellar ∼ 10^8 − 9 M_⊙, expected values according to the MBH–galaxy mass relation (see the left panel of Fig. 3). Interestingly, no redshift evolution is seen in the typical mass of the galaxy hosting LISA systems. In the same figure, we include the median stellar mass of galaxies harbouring single MBHs with the same mass as LISA binaries. This comparison enables us to determine whether or not galaxies housing LISA systems represent untypical hosts of low-mass MBHs. As shown, at fixed MBH mass, LISA systems are placed in slightly more massive galaxies than single MBHs. However, these differences are relatively small with values < 0.1 − 0.5 dex, depending on the specific mass of the binary. This trend is also seen in recent cosmological hydrodynamical simulations. For instance, the work of Dong-Páez et al. (2023a) showed that when post-processing dynamical delays between galaxies and MBH mergers are taken into account in the OBELISK simulation, merging MBHs tend to be placed slightly higher in the M_BH − M_stellar scaling relation than singles MBHs. Despite this, the difference is small enough that they agree within the scatter of the global relation.

Fig. 6. Median stellar mass of galaxies hosting LISA sources (red) and single MBHs (blue) at different redshifts. Red and blue areas represent the percentile 16th–84th. Different panels correspond to MBHs and MBHBs with masses: 104 − 105 M⊙ (top), 105 − 106 M⊙ (middle), and 106 − 107 M⊙ (bottom).

The results presented above suggest that LISA MBHBs will be hosted by dwarf galaxies. Given that these represent the most abundant population of galaxies in the Universe, it is fundamental to determine whether or not the LISA hosts display any distinctive property allowing for their unequivocal identification. Motivated by this, in Fig. 7 we explore the specific star formation rate (sSFR = SFR/M_stellar) of the LISA MBHBs hosts. As expected, their values decrease towards low z but they are systematically larger than 10⁻¹⁰ yr⁻¹, which is compatible with a population of star-forming galaxies. For comparison, Fig. 7 displays the sSFR evolution of galaxies hosting single MBHs with the same mass as LISA binaries. As shown, their values and trends are indistinguishable from the ones featured by LISA hosts (see similar results presented in Dong-Páez et al. 2023a, but for higher redshifts). Therefore, the stellar activity of dwarf galaxies will not be a good discriminant to unequivocally pinpoint the galaxies where LISA MBHBs reside.

Fig. 7. Properties of the galaxies hosting LISA MBHBs and single MBHs with the same mass as LISA sources. Left panel: median sSFR of galaxies hosting LISA sources (red) and single MBHs (blue) at different redshifts. The grey area represents the region where galaxies are considered passive galaxies (sSFR < 10−11 yr−1). Right panel: same as the left panel but for the gas fraction (fgas). The grey dotted line highlights the value fgas = 0.5. In all the panels, red and blue areas represent the 16th–84th percentile. The results shown in this figure highlight that the galaxies hosting LISA MBHBs are gas-rich and star-forming.

In addition to specific star formation, in the right panels of Fig. 7 we show the gas fraction of the LISA hosts, f_gas, defined as M_cold/(M_cold + M_stellar). As shown, LISA binaries are placed in galaxies with a large gas content (> 60%), regardless of redshift. Similar trends have been reported by Li et al. (2022), who, by analysing the outputs of the Illustris-TNG hydrodynamical simulation (Nelson et al. 2019), showed that most of the detected LISA MBHBs would be located in gas-rich galaxies with gas fractions in the range of 0.6 − 0.9.

In line with what was seen before, the hosts of 10⁴ − 10⁶ M_⊙ MBHBs do not show important redshift variations in f_gas. Conversely, the galaxies of 10⁶ − 10⁷ M_⊙ MBHBs tend to be less gas rich towards low z, which is probably a consequence of their larger stellar mass content (the growth of the stellar component implies a significant depletion of the gas reservoir). Following the comparison made with the stellar mass, Fig. 6 includes the gas fraction of galaxies harbouring single MBHs with the same mass as our selected LISA systems. Despite the absence of evident differences, the hosts of single MBHs display slightly larger f_gas. This is most likely due to their smaller stellar mass, suggesting that a smaller amount of their gas was converted into stars.

In addition to the intrinsic galaxy properties, such as star formation rate or gas fraction, in Fig. 8 we explore the possibility that galaxy morphology could be used as a guide when selecting low-mass galaxies harbouring LISA systems. Specifically, Fig. 8 shows the fraction of LISA MBHBs placed in ellipticals and spiral galaxies (divided into pseudobulges, classical bulges, and extreme late-type; see Sect. 2.3). As shown, at z > 1.5 ∼75% of MBHBs with M_Bin = 10⁴ − 10⁷ M_⊙ are hosted in extreme late-type systems. This trend is the result of the merger history of these high-z galaxies, which is ruled by interactions with small baryonic merger ratios (either extreme minor mergers or smooth accretion) whose capability of building bulges is negligible. Interestingly, at these high redshifts, elliptical and pseudobulge morphologies have a marginal contribution at any mass, representing less than 20% of the total population.

Fig. 8. Fractions of LISA MBHB sources (left panels) and single MBHs (right panels) hosted in elliptical (red), classical bulge (green), pseudobulge (blue), and extreme-late type (purple) galaxies at different redshifts. Top, middle, and bottom panels correspond to different mass bins: 104 − 105 M⊙, 105 − 106 M⊙, and 106 − 107 M⊙, respectively. We note that the MBHB hosts have noisier distributions because their number density is up to 2 − 3 orders of magnitude smaller than that of single MBHs (see Fig. 5). Briefly, galaxies hosting LISA systems display a disc-dominated morphology whose bulge component is more pronounced towards lower redshifts.

At z < 1.5, the trends for 10⁴ − 10⁷ M_⊙ and 10⁵ − 10⁶ M_⊙ change and approximately 75% of these MBHBs inhabit disc-dominated galaxies with a classical bulge. For systems with 10⁶ − 10⁷ M_⊙, classical bulges dominate as well, being the typical bulge structure for 50% of the population. However, pseudobulges and elliptical structures are also of considerable relevance, representing the other approximately 50% of the population. As a reference, in the left panels of Fig. 8 we show the morphological properties of the galaxies hosting single MBHs with the same mass as the LISA MBHBs. Interestingly, the morphology of these galaxies does not share the redshift evolution seen in the LISA hosts. At any redshift and mass, disc-dominated galaxies with an extreme late-type morphology dominate (40%−75%) the hosts of low-mass (< 10⁷ M_⊙) single MBHs. This small difference seen between the morphology of normal dwarf galaxies and LISA hosts could help in the identification of the galaxies where z < 1 MBHBs of 10⁴ − 10⁶ M_⊙ are placed: low-mass galaxies with a more predominant bulge component than the average population are more likely to harbour a LISA system.

Finally, we explore the optical counterpart of the galaxies where LISA sources reside. Our aim is to determine the possibility of detecting the LISA host in cases where no AGN counterpart associated with the MBHB is found, that is, in cases where the binary is in an inactive phase. To this end, in Fig. 9 we present the apparent magnitude of the LISA hosts in the optical r band³. We stress that similar behaviours are found in the other optical bands such as g or i. As expected, high-z galaxies are dimmer than their low-z counterparts. For instance, z > 2 galaxies have r ∼ 27, while galaxies at z < 0.5 display r < 25. In addition, Fig. 9 shows that the magnitude of LISA MBHB hosts does not increase (i.e become dimmer) at z > 1 but remains constant. This is the effect of the fast rise of the galaxy star formation towards high z (see Fig. 7), which is able to compensate for the effect of the luminosity distance in making sources dimmer. In the same figure, we include the median r band magnitude of galaxies harbouring single MBHs with the same mass as LISA binaries. As shown, no significant differences are seen between these two samples except for the case of 10⁶ − 10⁷ M_⊙, in which the hosts of LISA systems are slightly brighter than the hosts of single MBHs. This deviation is driven by the stellar mass of the galaxy hosts, which tend to be slightly larger for the case of MBHBs (see Fig. 6).

Fig. 9. Observability of the galaxies hosting LISA MBHBs. Left panel: median r-band magnitude of the galaxies hosting LISA MBHBs (red) and single MBHs (blue). Shaded areas represent the percentile 16th–84th. Dashed grey lines correspond to the r band limiting magnitude of SDSS (r = 25.1) and LSST (r = 27.5). Right panel: fraction of galaxies hosting LISA MBHBs (red) and single MBHs (blue) whose r-band magnitude is smaller than rth = 25.1 (solid line) and 27.5 (dashed line), i.e. those that are detectable. Top, middle, and bottom panels correspond to different mass bins: 104 − 105 M⊙, 105 − 106 M⊙, and 106 − 107 M⊙, respectively. Overall, the figure shows that LISA MBHBs are hosted in dim galaxies; indeed, half of them are fainter than the detection limit of current photometric surveys.

To guide the reader as to the observability of LISA MBHB host galaxies, in Fig. 9 we highlight the r band detection limits of SDSS (r = 25.1) and LSST (r = 27.5). As shown, SDSS will only be able to detect the optical emission of z < 0.5 LISA hosts. On the other hand, LSST can extend this detection up to z ∼ 2 − 2.5. Consequently, these results highlight that without any AGN emission raised by the MBHB, the optical identification of LISA hosts will only be feasible at z ≤ 1, while at higher redshift the galaxies are at the limits of detection capabilities. To show this, the redshift evolution of the fraction of galaxies below the detection limit of SDSS and LSST (i.e not detectable) is presented in the right panel of Fig. 9. As shown, at z > 1, only around 20% and 50% of the galaxies will be found by the optical surveys SDSS and LSST, respectively.

Overall, the results presented in this section indicate that LISA MBHBs will be placed in dim galaxies with small stellar content, large gas fraction, and an active star formation history. However, they will not be peculiar systems, because few differences are found with respect to the galaxies housing single MBHs with the same mass as LISA binaries. In this way, the unequivocal identification of LISA hosts through standard galaxy properties will be challenging. Motivated by this, in the following section, we explore the presence of merger signatures as a tracer that can be used to select LISA hosts among the population of dwarf galaxies.

5. Merger signatures: An important feature to detect the hosts of LISA MBHBs?

According to our current paradigm, galaxy mergers are essential requisites for the creation of MBHBs (Begelman et al. 1980). In accordance, galaxies hosting LISA systems might display visible merger signatures of the galaxy interaction that led to the MBHB formation. In this section, we explore the feasibility of pinpointing LISA MBHB hosts through the identification of merger signatures in dwarf galaxies. We stress that the merger ratios reported in this work correspond to the baryonic ones, that is, they account for the stellar and gas component.

5.1. Merger signatures in the hosts of LISA MBHBs

The left panel of Fig. 10 presents the baryonic mass ratio of the galaxy interaction that deposited the secondary MBH to the galaxy, hereafter $m_{\mathrm{R}}^{\mathrm{G},\sec }$. For simplicity, we only present the LISA MBHBs at three different redshifts bins: z = 1 ± 0.25, 2 ± 0.25, and 3 ± 0.25. As shown, 50% of the mergers that brought the secondary MBH to the galaxy display $0.002<m_{\mathrm{R}}^{\mathrm{G},\sec }<0.02$, regardless of redshift and binary mass. On the other hand, major mergers ($m_{\mathrm{R}}^{\mathrm{G},\sec }>0.2$) contribute 10% of the cases.

Fig. 10. Merger history of the hosts of LISA MBHBs. Left panel: cumulative distribution of the baryonic mass ratio ($m_{\mathrm{R}}^{\mathrm{G},\sec }$) of the two galaxies involved in the interaction that brought the secondary MBH to the primary galaxy. Different columns correspond to different redshift bins: z = 1 ± 0.25 (left), z = 2 ± 0.25 (middle), and z = 3 ± 0.25 (right). The horizontal dashed grey lines highlight the values 0.1, 0.5, and 0.9. The shaded region corresponds to the merger ratios associated with major mergers in L-Galaxies. Right panel: same as left panel but for the values of the delay time, ΔtBin, between the galaxy–galaxy merger and that of the binary. In this case, the shaded area represents the region where the merger signatures have vanished at the time of LISA detection. Dashed and dotted lines represent the same but for galaxy interactions with mR > 0.2 and 0.01, respectively. In all panels, vertical dashed lines correspond to the Hubble time at z = 1 (left panel), z = 2 (middle panel), and z = 3 (right panel). In general, the distributions shown in these figures show that more than 50% of the LISA MBHBs result from galaxy mergers with a baryonic mass ratio of greater than 0.01, with a minor contribution to major mergers. We stress that these baryonic ratios do not imply that the binaries display the same mass ratios at inspiral and coalescence. This is due to the gas accretion of MBHBs along the galaxy evolution.

Besides merger ratios, another important quantity to take into account is the moment at which these mergers took place. This quantity can provide precious information about still visible signatures related to the galaxy interaction, such as stellar tidal tails, bridges, streams, and shell structures (Toomre & Toomre 1972; Gerber & Lamb 1994; Lotz et al. 2004). For instance, by analysing hydrodynamical cosmological simulations, Mancillas et al. (2019) found that major and intermediate merger events (i.e. stellar merger ratios ≥0.1) leave post-merger signatures with a survival timescale of 0.7 − 4 Gyr. While tidal tails remain visible between 0.7 and 1 Gyr, shells and streams could last up to 3 − 4 Gyr and 1.5 − 3 Gyr, respectively. Taking into account these results, hereafter we assume that merger signatures last on average 2 Gyr, independently of the high redshift explored here. We highlight that this duration of merger signatures should be considered as an upper limit for high-redshift galaxies. As the galaxy dynamical time goes as (1 + z)^−3/2 (Mo et al. 1998), the time period over which distortion or non-asymmetric features brought about by mergers will prevail will be up to eight times shorter at z = 3 than at z = 0 (see further discussion in Volonteri et al. 2021). On top of this, we stress that we do not take into account the surface brightness of these features and the results presented in this section should be considered as upper limits; that is, the merger signature(s) could be present in the galaxy but its associated surface brightness could be sufficiently low to hamper its detection.

To explore the presence of merger signatures related to the galaxy interaction that brought the secondary MBH to the galaxy, in Fig. 10 we present the cumulative distribution of Δt_Bin, defined as:

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{Bin}}=t_{\mathrm{now}}-t_{\mathrm{G},\sec },\end{array}$$(7)

where t_now is the lookback time associated with the redshift at which the binary is detected, and t_G, sec is the lookback time at which the galaxy merger leading to the binary formation took place. To guide the reader, the smaller the value of Δt_Bin, the smaller the time elapsed between the observation (i.e. GW detection) and the galaxy merger involved in the formation of the binary. As shown, the large majority of LISA hosts at z ∼ 3 display Δt_Bin < 1 Gyr, suggesting that they would still display merger signatures at the moment of the MBHB detection. Galaxies at z ∼ 3 hosting 10⁶ − 10⁷ M_⊙ are the ones that have the smallest Δt_Bin values, with 50% of them having Δt_Bin < 0.8 Gyr. These small Δt_Bin values at z ∼ 3 are caused by the fact that MBHB evolution (pairing and hardening) is faster for dense high-z galaxies with relatively massive MBHs.

When the distributions are divided by merger ratios, the galaxy interactions associated with major mergers (m_R > 0.2) are skewed towards smaller Δt_Bin values irrespective of the mass bin studied. This shows that, as expected, pairing and hardening evolution is faster in major galaxy mergers. Despite the fact that some displacement can also be seen for intermediate-merger ratios (m_R > 0.01) in M_Bin > 10⁵ − 10⁷ M_⊙, these are systematically smaller than in the case of major interactions. According to these distributions, 60%−70% of the galaxies hosting z ∼ 3 LISA sources would display a merger signature related to the major or intermediate interaction that led to the MBHB formation. These described trends are kept when lower redshifts are studied. However, fewer galaxies housing LISA sources would display merger signatures. For instance, at z ∼ 2 (z ∼ 1) only 50% (30%−40%) of the LISA sources associated with major or intermediate galaxy mergers display signs that reveal the interaction, regardless of the MBHB mass. In particular, among all the galaxies harbouring LISA systems, the ones with M_Bin = 10⁴ − 10⁵ M_⊙ at z ∼ 1 are the ones that show the least presence of merger signatures related to the formation of the MBHB, with < 10% of cases showing such signatures.

Our analysis of the assembly history of LISA hosts highlights the fact that at the moment of the GW detection, some galaxies harbouring LISA sources would still display merger signatures related to the interaction that brought the secondary MBH to the galaxy. Despite this, the lifetime of a galaxy is rather complicated and multiple mergers can happen between the observation time and the merger responsible for the binary formation (see Volonteri et al. 2020). To illustrate this, in Fig. 11 we present the merger tree of two random galaxies hosting 10⁵ − 10⁶ M_⊙ and 10⁶ − 10⁷ M_⊙ LISA MBHBs at z ∼ 1. As we can see, the galaxy interaction responsible for the MBHB formation took place at z ∼ 12 − 11 but the galaxy underwent subsequent galaxy mergers until z ∼ 1. Therefore, the LISA hosts might show merger signs not correlated with the ones raised by the galaxy interaction leading to the MBHB formation. This can cause confusion in cases where the merger features are used to infer the mass of the satellite galaxy involved in the interaction or the time at which that merger took place and therefore determine the binary lifetime.

Fig. 11. Galaxy merger tree of the LISA host at z ∼ 1.0. The time evolution is represented as a function of the lookback time (tlookback). For reference, tlookback = 8.5 (11.5) Gyr corresponds to z ∼ 1 (z ∼ 13). Vertical black lines connect the secondary branches with the main ones (always at the bottom). Thus, the linking points between branches correspond to galaxy mergers. We highlight the baryonic mass ratio of the merging galaxies that caused the binary formation in red. The colour of each dot encodes the stellar mass of the galaxy. Finally, the extension of the grey arrow corresponds to the length of ΔtBin. The two galaxy merger trees reported in this figure show that multiple episodes of galaxy mergers can happen before and after the formation of the LISA source.

To explore the level of confusion, in Fig. 12 we present the probability, 𝒫, that a LISA MBHB host experienced a further merger within Δt_Bin. As shown, regardless of the MBHB mass, LISA hosts have 𝒫 values of > 80%. This probability drops when different galaxy merger ratios are taken into account. Concerning major mergers (m_R > 0.2), the hosts of 10⁴ − 10⁷ M_⊙ MBHBs placed at z > 1.5 display a ∼10% probability of undergoing one of these events within Δt_Bin. This value rises up to 40%−50% when lower redshifts are considered.

Fig. 12. Probability, 𝒫, that a galaxy hosting a LISA MBHB (left) and a single MBHs (right) will undergo a galaxy merger within ΔtBin. Solid, dashed, and dotted lines correspond to any baryonic mass ratio (mR), mR > 0.2, and mR > 0.01, respectively. Top, middle, and bottom panels: results for MBin = 104 − 105 M⊙, 105 − 106 M⊙, and 106 − 107 M⊙, respectively. The large probability values shown in this figure show that galaxies hosting the LISA source display merger signatures that are uncorrelated with the galaxy merger that led to the MBHB formation and coalescence.

When accounting for intermediate merger ratios (m_R > 0.01), the results show a similar trend to the one seen in the case of major mergers. However, some differences are seen for the hosts of 10⁴ − 10⁵ M_⊙ and 10⁵ − 10⁶ M_⊙ MBHBs where 𝒫 can be a factor of about 3 larger at high z, reaching values of 40%. Finally, in the right panel of Fig. 12 we present the probability 𝒫 for systems harbouring single MBHs with the same mass as our selected LISA MBHBs. In these cases, it is not possible to set the value of t_G, sec (see Eq. (7)). As a reference, we have chosen the median t_G, sec computed from all the LISA MBHBs placed at the explored redshift. The figure shows small differences for major merger events, where single MBHs at z < 1 display values of 𝒫 that are around 10% smaller. Nevertheless, the trends and values of 𝒫 are shared between the LISA and single MBH hosts.

The analysis performed in this section of the merger history of LISA hosts reveals that merger features should be taken with caution. The galaxies where LISA MBHBs are placed will often display merger signatures that are uncorrelated with the galaxy interaction that led to the LISA MBHB formation. In addition, dwarf galaxies not hosting MBHBs can also display similar merger features.

5.2. Merger signatures in the fields of LISA MBHBs

Recent works have shown that once the GW signal is detected, LISA can localise the MBHB with a sky area that ranges from several hundreds of deg² to fractions of a deg², depending on the intrinsic properties of the binary and its redshift. Given these potentially wide areas, the number of candidates can be relatively large, hampering the unequivocal detection of the LISA host galaxy. Indeed, using simulated data, Lops et al. (2023) showed that the number of galaxies lying within the LISA sky area can be as large as 10⁵. This implies that the unequivocal detection of LISA hosts would benefit from a pre-selection based on some specific galaxy properties or features. Driven by this, in this section, we explore how feasible it is to identify the LISA host among all the galaxies in its surroundings using merger features. Hereafter, we define the LISA MBHB environment as the distribution of galaxies within 8 Mpc. This distance selection has been chosen so that the sky-projected area resembles the sky localisation of LISA when it detects a z > 1 MBHB at merger time (∼0.1 − 0.01 deg²; see e.g. Mangiagli et al. 2020; Piro et al. 2023). In addition to this definition, we focus on galaxies with a stellar mass of 10⁸ < M_stellar < 10⁹ M_⊙, that is, systems with stellar masses similar to those of LISA hosts (see Sect. 4).

Figure 13 presents the environment of two random LISA binaries with masses of 10⁵ − 10⁶ M_⊙ and 10⁶ − 10⁷ M_⊙ at z = 2. As we can see, the number of galaxies with similar stellar mass to the LISA hosts is relatively large. When these galaxies are highlighted according to their merger history, a large number of them experienced a galaxy interaction with an intermediate merger ratio of (0.01 < m_R < 0.2) within the last 2 Gyr (i.e. the average lifetime of merger features). On the other hand, the number of galaxies that underwent a major merger (m_R > 0.2) within the last 2 Gyr is smaller but not negligible. These plots suggest that the environment of LISA sources would be crowded by dwarf galaxies featuring merger signatures. This would hinder the unequivocal detection of LISA sources by pre-selecting galaxies with recent signs of interaction.

Fig. 13. Sky map of a LISA source at z ∼ 2 in the simulated universe. The left panel displays the case of a MBin = 105 − 106 M⊙ LISA MBHB, and the right panel corresponds to a MBHB of 106 − 107 M⊙. Each dot corresponds to a galaxy with a mass of 108 < Mstellar < 109 M⊙ whose distance to the LISA source is smaller than 8 Mpc. Black dots represent galaxies that never experienced a merger or a merger that happened more than 2 Gyr ago. Red dots correspond to galaxies that experienced at least one major merger in the last 2 Gyr. Blue dots display the galaxies that experience a merger with a mass ratio of 0.01 < mR < 0.2 in the last 2 Gyr.

To quantify the previous results, in Fig. 14 we present the probability 𝒯 that a galaxy with 10⁸ < M_stellar < 10⁹ M_⊙ within 8 Mpc of the LISA host displays a merger signature⁴. We explore three different redshift bins: z = 1 ± 0.25, 2 ± 0.25, and 3 ± 0.25. As shown, the probability 𝒯 displays a redshift evolution regardless of the considered mass bin. At z ∼ 3, the value of 𝒯 reaches 100% while at z ∼ 1 the probability drops down to about 80%. This redshift dependence vanishes when only moderate baryonic merger ratios are considered (i.e. m_R > 0.01). Specifically, a rather constant probability of about 60% is reached in these cases. Finally, for dwarf galaxies displaying major merger (m_R > 0.2) signatures, the results show an increasing probability towards low z. Whereas at z ∼ 3 galaxies with major merger signatures have an approximately 10% chance of appearing in the LISA host environments, at z ∼ 1 the value of 𝒯 increases up to roughly 25%.

Fig. 14. Probability, 𝒯, that the environments of z = 1 ± 0.25, 2 ± 0.25, and 3 ± 0.25 LISA MBHBs have at least one galaxy of 108 < Mstellar < 109 M⊙ featuring merger signatures. Black, red, and blue dots correspond to the values of 𝒯 when the population of 108 < Mstellar < 109 M⊙ galaxies is divided according to galaxy merger ratios: all ratios, mR > 0.2 (major mergers) and mR > 0.01 (intermediate mergers), respectively. The error bars correspond to the 32th–68th percentile. Horizontal dashed lines highlight the values 𝒯 = 0.1, 0.5, and 0.9.

6. Caveats

In this section, we discuss some caveats related to the model scheme of MBHs and MBHBs that can cause variations in the results presented in this work. As shown in Izquierdo-Villalba et al. (2020), the physical prescriptions of L-Galaxies applied on top of the Millennium merger trees generate a population of active MBHs that resembles the one reported in observations. However, Spinoso et al. (2023) showed that L-Galaxies tend to over-predict the growth of high-z (z > 3) small MBHs when the higher-resolution merger trees of Millennium-II are used and the MBH spin is neglected. In fact, this shortcoming is not a unique feature of L-Galaxies, and other semi-analytical models and hydrodynamical simulations feature it (see e.g. Degraf et al. 2010; DeGraf & Sijacki 2020; Marshall et al. 2020; Trinca et al. 2022). This over-prediction could cause the model to more frequently generate over-massive black holes in low-mass galaxies. Consequently, the predictions concerning the stellar content of LISA hosts should be considered as lower limits. In a future paper (Spinoso et al., in prep.), we will perform a further investigation of MBH growth in dwarf galaxies.

Another caveat to take into account is the specific treatment of the hardening phase of MBHBs. Since L-Galaxies cannot track the stellar distribution inside galactic bulges, assumptions of stellar profiles are required. Specifically, our SAM uses the Sérsic model, which has been shown to be a good profile to represent the bulge population at several cosmological times (see e.g. Drory & Fisher 2007; Fisher & Drory 2008; Gadotti 2009; Shibuya et al. 2015). However, as proved by Biava et al. (2019), different assumptions about the stellar mass distributions around the binaries lead to variations in the MBHB lifetimes. For instance, the commonly used Dehnen profile (Dehnen 1993) predicts a smaller MBHB lifetime than the Sérsic model used in this work. On a related topic, the L-Galaxies model neglects the presence of nuclear stellar clusters (a feature that will be addressed in Polkas et al., in prep.; Hoyer et al., in prep.). Including the large stellar concentrations of this kind of system would lead to a shorter hardening phase of the MBHBs, changing the expected distribution of the number density of LISA systems predicted in this work.

7. Conclusions

In this paper, we study the properties of galaxies hosting LISA MBHBs with masses of 10⁴ − 10⁷ M_⊙ at z ≤ 3. To this end, we generated a simulated lightcone using the L-Galaxies semi-analytical model applied to the high-resolution Millennium-II merger trees (5.7 × 10⁷ − 3 × 10¹⁴ M_⊙ halo mass range). The version of the SAM used in this work corresponds to that presented in Izquierdo-Villalba et al. (2020, 2022a) with an improved prescription for the formation of MBHs (based on Spinoso et al. 2023). This L-Galaxies variant includes different physical models in order to self-consistently tackle the growth of MBHs and the dynamical evolution of MBHBs (from the galaxy merger down to the GW inspiral phase). The resulting lightcone contains galaxies, MBHs, and MBHBs of different masses up to z ∼ 3, and features a line of sight (RA, Dec) = (77.1, 60.95) deg with an angular extension of (δRA, δDec) = (45.6, 22.5) deg.

The main results of this work can be summarised as follows:

• LISA systems represent less than 1% of the 10⁴ − 10⁷ M_⊙ MBH population, and their abundance is not strongly correlated with the mass of the binary. While at z ∼ 3, LISA systems with a total mass of 10⁴ − 10⁷ M_⊙ display abundances of ∼0.1 deg⁻², at z ∼ 0.5 they have ∼0.01 deg⁻².

• LISA MBHBs are hosted in dwarf galaxies with stellar masses of 10⁸ − 10⁹ M_⊙ at any z < 3. At fixed total black hole mass, these hosts are ∼0.1 − 0.5 dex more massive than the ones housing single MBHs. However, these differences are not large enough to be easily measured from observations.

• The galaxies hosting LISA MBHBs are gas-rich (gas fraction > 0.6) and star-forming galaxies (sSFR > 10⁻¹⁰ yr⁻¹). Galaxies hosting single black holes of comparable mass share similar properties.

• LISA hosts have a disc-dominated morphology whose spheroidal component varies with redshift. While systems at z > 1 display a tiny or negligible bulge contribution (B/T < 0.01), at lower redshifts the spheroidal component plays a more important role (0.01 < B/T < 0.7). Galaxies harbouring single MBHs with the same mass as LISA sources are also placed in disc-dominated galaxies. However, their bulge importance does not show any redshift evolution, showing a negligible contribution to the galaxy morphology (B/T < 0.01).

• LISA MBHBs are placed in faint galaxies, with a magnitude in the r optical band of greater than 20. Taking this into account, as well as the fact that the MBHB might be in an inactive phase (i.e no AGN emission), current optical facilities will only be able to identify the LISA hosts at z ≤ 0.5. The use of surveys like LSST will extend this search up to z ∼ 2.

The small differences found (at a fixed MBH mass) between the hosts of LISA MBHBs and single MBHs underline that the unequivocal identification of LISA hosts through standard galaxy properties will be challenging. Motivated by this, we also explored the presence of merger signatures as a property that can be used to select LISA hosts among the population of dwarf galaxies. Our main results on this topic can be summarised as follows:

• Around 80% of the secondary MBHs forming LISA MBHBs were deposited in the galaxy after the accretion of a companion galaxy with approximately 5 − 100 times smaller baryonic mass. These interactions left host merger signatures in the LISA galaxy, which are visible for 80% of the cases at z ∼ 2. This fraction drops down to 40% for z ∼ 1 hosts.

• Given the long lifetime of LISA MBHBs, their hosts have a high chance of undergoing several galaxy interactions. This implies that LISA hosts will have a large probability of displaying merger signatures uncorrelated with the ones produced by the merger that led to the binary formation. In particular, we find that this probability mildly depends on the binary mass and rises towards low z, with values reaching about 60%−80% at z < 1 when accounting for mergers with a baryonic mass ratio of > 0.01.

• The environments of LISA hosts will be populated by a large number of dwarf galaxies displaying signs of recent mergers as well. Indeed, about 60% of them will have merger signatures caused by galaxy interactions with an intermediate baryonic merger ratio (> 0.01), that is, similar to the ones displaying LISA hosts.

Taking into account all the results summarised above, LISA hosts can be expected to be faint, low-mass galaxies whose intrinsic characteristics – such as star-formation activity – or extrinsic properties – such as merger signatures – will not be good tracers for their unequivocal identification. In view of these results, the proposal and study of alternative approaches for pinpointing LISA hosts would be welcome. For instance, strategies based on the search for EM counterparts associated with faint AGNs with light curves and spectra characteristic of a merger (d’Ascoli et al. 2018; Yuan et al. 2021) may be promising in this respect (Mangiagli et al. 2022; Lops et al. 2023; Dong-Páez et al. 2023b).

Acknowledgments

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. M.V. acknowledges funding from the French National Research Agency (grant ANR-21-CE31-0026, project MBH_waves) and from the Centre National d’Études Spatiales. 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 partial support from the project PGC2018-097585-B-C22.

References

All Figures

	Fig. 1. Probability 𝒫 that a given DM halo of mass Mhalo would be seeded by a MBH. Each color corresponds to a different redshift. Vertical and horizontal grey dashed lines highlight the halo mass resolution of MSII and the value 𝒫 = 0.5, respectively.
In the text

	Fig. 2. Massive black hole occupation fraction at z = 0 as a function of stellar mass (Mstellar). The shaded grey region corresponds to the constraints presented in the observational work of Miller et al. (2015).
In the text

Fig. 3. General properties of MBHs in L-Galaxies SAM applied in the Millennium-II merger trees. Left panel: 2D histogram of the bulge-black hole correlation at z = 0. The predictions are compared with the observations of Kormendy & Ho (2013; K&H2013) where blue, orange, and red points correspond to pseudobulges, classical bulges, and ellipticals, respectively. Right panel: black hole mass function at z = 0 when the model is run on top of the Millennium-II merger trees (black lines). Error bars correspond to the Poissonian error. For comparison, we have added the observational constraints of Marconi et al. (2004) and Shankar et al. (2004, 2009, 2013).

In the text

Fig. 4. Quasar bolometric luminosity functions (Lbol) at z ∼ 0.5, 1.0, 2.0, and 3.0 predicted by L-Galaxies when MS (red) and MSII (black) DM merger trees are used. For comparison, the grey dashed lines represent the luminosity functions of L-Galaxies applied in addition to MSII when a cut of ≥107 M⊙ in black hole mass is performed. For comparison, these predictions are compared with the data presented in Shen et al. (2020; blue circles). In all the plots, error bars correspond to the Poissonian error (large error bars correspond to the cases where only one object is found).

In the text

Fig. 5. Number of hard MBHBs (grey solid) and single MBHs (grey dashed) per deg2 (dN(z)/dΩ, black) as a function of redshift. Different colors represent the same but when the population of MBHBs is divided by semi-major axis (aBH). Each panel corresponds to binaries with different total masses: 104 − 105 M⊙ (top), 105 − 106 M⊙ (middle), and 106 − 107 M⊙ (bottom). Overall, the figure shows that LISA MBHBs represent 0.1%−1% of the MBHs between 104 and 107 M⊙.

In the text

	Fig. 6. Median stellar mass of galaxies hosting LISA sources (red) and single MBHs (blue) at different redshifts. Red and blue areas represent the percentile 16th–84th. Different panels correspond to MBHs and MBHBs with masses: 104 − 105 M⊙ (top), 105 − 106 M⊙ (middle), and 106 − 107 M⊙ (bottom).
In the text

Fig. 7. Properties of the galaxies hosting LISA MBHBs and single MBHs with the same mass as LISA sources. Left panel: median sSFR of galaxies hosting LISA sources (red) and single MBHs (blue) at different redshifts. The grey area represents the region where galaxies are considered passive galaxies (sSFR < 10−11 yr−1). Right panel: same as the left panel but for the gas fraction (fgas). The grey dotted line highlights the value fgas = 0.5. In all the panels, red and blue areas represent the 16th–84th percentile. The results shown in this figure highlight that the galaxies hosting LISA MBHBs are gas-rich and star-forming.

In the text

Fig. 8. Fractions of LISA MBHB sources (left panels) and single MBHs (right panels) hosted in elliptical (red), classical bulge (green), pseudobulge (blue), and extreme-late type (purple) galaxies at different redshifts. Top, middle, and bottom panels correspond to different mass bins: 104 − 105 M⊙, 105 − 106 M⊙, and 106 − 107 M⊙, respectively. We note that the MBHB hosts have noisier distributions because their number density is up to 2 − 3 orders of magnitude smaller than that of single MBHs (see Fig. 5). Briefly, galaxies hosting LISA systems display a disc-dominated morphology whose bulge component is more pronounced towards lower redshifts.

In the text

Fig. 9. Observability of the galaxies hosting LISA MBHBs. Left panel: median r-band magnitude of the galaxies hosting LISA MBHBs (red) and single MBHs (blue). Shaded areas represent the percentile 16th–84th. Dashed grey lines correspond to the r band limiting magnitude of SDSS (r = 25.1) and LSST (r = 27.5). Right panel: fraction of galaxies hosting LISA MBHBs (red) and single MBHs (blue) whose r-band magnitude is smaller than rth = 25.1 (solid line) and 27.5 (dashed line), i.e. those that are detectable. Top, middle, and bottom panels correspond to different mass bins: 104 − 105 M⊙, 105 − 106 M⊙, and 106 − 107 M⊙, respectively. Overall, the figure shows that LISA MBHBs are hosted in dim galaxies; indeed, half of them are fainter than the detection limit of current photometric surveys.

In the text

Fig. 10. Merger history of the hosts of LISA MBHBs. Left panel: cumulative distribution of the baryonic mass ratio ($m_{\mathrm{R}}^{\mathrm{G},\sec }$) of the two galaxies involved in the interaction that brought the secondary MBH to the primary galaxy. Different columns correspond to different redshift bins: z = 1 ± 0.25 (left), z = 2 ± 0.25 (middle), and z = 3 ± 0.25 (right). The horizontal dashed grey lines highlight the values 0.1, 0.5, and 0.9. The shaded region corresponds to the merger ratios associated with major mergers in L-Galaxies. Right panel: same as left panel but for the values of the delay time, ΔtBin, between the galaxy–galaxy merger and that of the binary. In this case, the shaded area represents the region where the merger signatures have vanished at the time of LISA detection. Dashed and dotted lines represent the same but for galaxy interactions with mR > 0.2 and 0.01, respectively. In all panels, vertical dashed lines correspond to the Hubble time at z = 1 (left panel), z = 2 (middle panel), and z = 3 (right panel). In general, the distributions shown in these figures show that more than 50% of the LISA MBHBs result from galaxy mergers with a baryonic mass ratio of greater than 0.01, with a minor contribution to major mergers. We stress that these baryonic ratios do not imply that the binaries display the same mass ratios at inspiral and coalescence. This is due to the gas accretion of MBHBs along the galaxy evolution.

In the text

Fig. 11. Galaxy merger tree of the LISA host at z ∼ 1.0. The time evolution is represented as a function of the lookback time (tlookback). For reference, tlookback = 8.5 (11.5) Gyr corresponds to z ∼ 1 (z ∼ 13). Vertical black lines connect the secondary branches with the main ones (always at the bottom). Thus, the linking points between branches correspond to galaxy mergers. We highlight the baryonic mass ratio of the merging galaxies that caused the binary formation in red. The colour of each dot encodes the stellar mass of the galaxy. Finally, the extension of the grey arrow corresponds to the length of ΔtBin. The two galaxy merger trees reported in this figure show that multiple episodes of galaxy mergers can happen before and after the formation of the LISA source.

In the text

Fig. 12. Probability, 𝒫, that a galaxy hosting a LISA MBHB (left) and a single MBHs (right) will undergo a galaxy merger within ΔtBin. Solid, dashed, and dotted lines correspond to any baryonic mass ratio (mR), mR > 0.2, and mR > 0.01, respectively. Top, middle, and bottom panels: results for MBin = 104 − 105 M⊙, 105 − 106 M⊙, and 106 − 107 M⊙, respectively. The large probability values shown in this figure show that galaxies hosting the LISA source display merger signatures that are uncorrelated with the galaxy merger that led to the MBHB formation and coalescence.

In the text

Fig. 13. Sky map of a LISA source at z ∼ 2 in the simulated universe. The left panel displays the case of a MBin = 105 − 106 M⊙ LISA MBHB, and the right panel corresponds to a MBHB of 106 − 107 M⊙. Each dot corresponds to a galaxy with a mass of 108 < Mstellar < 109 M⊙ whose distance to the LISA source is smaller than 8 Mpc. Black dots represent galaxies that never experienced a merger or a merger that happened more than 2 Gyr ago. Red dots correspond to galaxies that experienced at least one major merger in the last 2 Gyr. Blue dots display the galaxies that experience a merger with a mass ratio of 0.01 < mR < 0.2 in the last 2 Gyr.

In the text

Fig. 14. Probability, 𝒯, that the environments of z = 1 ± 0.25, 2 ± 0.25, and 3 ± 0.25 LISA MBHBs have at least one galaxy of 108 < Mstellar < 109 M⊙ featuring merger signatures. Black, red, and blue dots correspond to the values of 𝒯 when the population of 108 < Mstellar < 109 M⊙ galaxies is divided according to galaxy merger ratios: all ratios, mR > 0.2 (major mergers) and mR > 0.01 (intermediate mergers), respectively. The error bars correspond to the 32th–68th percentile. Horizontal dashed lines highlight the values 𝒯 = 0.1, 0.5, and 0.9.

In the text