© ESO 2022

1. Introduction

Massive black holes (MBHs) residing in the centers of many (if not all) galaxies are fascinating objects which elude our basic understanding even after decades of study. Indeed, numerous observational campaigns have contributed strong evidence to support their existence in the nuclei of both large galaxies (Ferrarese & Ford 2005) and also dwarf galaxies (Reines et al. 2013; Nguyen et al. 2019; Baldassare et al. 2020), while also showing complex empirical relations between MBH masses and various properties of their host galaxies such as velocity dispersion, total stellar mass, and luminosity (Gültekin et al. 2009; Beifiori et al. 2012; McConnell & Ma 2013; Läsker et al. 2014; Reines & Volonteri 2015). An exact explanation of this apparently tight coevolution is still missing, although several hypotheses have been proposed. For instance, there might be one specific factor that has a strong influence on both galaxy and its MBH, such as the total gas reservoir (see Kulier et al. 2015 and references therein). Another compelling idea is that there is a mutual feedback: galaxy evolution processes regulate the inflow of the material towards its center (thereby the growth of the MBH) and, on the other hand, a rapidly accreting MBH can lead to gas heating and outflows, resulting in suppressed star formation (Maiolino et al. 2012).

One of the fundamental problems in the field is related to the masses of MBHs found in the range of 10⁴ − 10¹⁰ M_⊙ with the most massive and luminous quasars already being observed at redshifts as high as z = 7.54 (Wu et al. 2015; Mazzucchelli et al. 2017; Bañados et al. 2018). Firstly, it is still unclear what type of objects and processes constitute the formation of MBH seeds. There are currently several proposed scenarios, three of which have gained particular interest and are being extensively studied in recent years (Volonteri et al. 2021). These include (i) low mass seeds (∼100 M_⊙) from population III stars (PopIII), (ii) intermediate mass seeds (∼10⁴ M_⊙) from runaway collisions in stellar clusters and (iii) heavy seeds (10⁴ − 10⁵ M_⊙) from direct collapse black holes (for an extensive review see Inayoshi et al. 2020). The second unknown corresponds to processes responsible for the mass growth and feedback exerted back on the host galaxy. In this case, the main three channels are (Croton et al. 2006): (i) quasar (radiative) mode characterized by high accretion rates and strong X-ray emission (induced by galaxy mergers and/or disc instabilities), (ii) radio (jet) mode where accretion rates are lower due energy transfer to a jet production and (iii) MBH mergers.

It is extremely challenging to test the above models with electromagnetic (EM) observations. This is due to the fact that MBH evolutionary processes started very early in the history of the Universe and the most distant objects that we can now study represent only a small and highly biased fraction of the whole population (Barausse et al. 2017; Shankar et al. 2019). Theoretical predictions indicate that MBH seeds from remnants of PopIII stars could form at redshifts as high as 20 (Volonteri 2010; Inayoshi et al. 2020), while redshifts considered for DCBHs fall approximately between 13 and 20 (Yue et al. 2014). This is outside the range of any of current EM observatories, however future missions like SKA in radio (Taylor 2000), JWST in IR/optical (Behroozi et al. 2020), and Athena in X-rays (Nandra et al. 2013) will probe sources up to redshifts of about 15−20.

Both observations and theory predict that MBHs can form binaries and merge within the Hubble time. It is one of the fundamental assumptions in hierarchical models of galaxy evolution (Benson 2010) and there is a non-negligible number (a few hundred) of observed candidate systems that show signs of hosting an inspiraling MBH binary at pc–subpc separations (MBHB; De Rosa et al. 2019; Saade et al. 2020). Such systems are expected to emit low-frequency gravitational waves (GWs) in the range of about nHz to mHz, depending on binary mass and separation. As GWs hardly interact with matter, they can carry unperturbed information further than EM radiation and thus may enable us to reach the earliest stages of the Universe’s evolution.

LISA (Laser Interferometer Space Antenna, Amaro-Seoane et al. 2017) is scheduled to start its operation around 2034 and will be sensitive to mergers of MBHBs of around 10³ − 10⁷ M_⊙, radiating at mHz frequencies up to redshifts well above 20. Moreover, lower (nHz) part of the GW spectrum is the main target for Pulsar Timing Arrays (PTAs) already being in operation: PPTA (Manchester et al. 2013), EPTA (Desvignes et al. 2016), and NANOGrav (McLaughlin 2013), as well as (most recently) inPTA (Joshi et al. 2018) and CPTA (Lee 2016). What is particularly important is that PTAs have also already constrained an upper limit for the GW background (GWB) amplitude generated by an assembly of cosmological population of MBHBs to approximately 10⁻¹⁴ <  A <  10⁻¹⁵ at a reference frequency of f = 1/yr (Shannon et al. 2015; Babak et al. 2016; Arzoumanian et al. 2018). These measurements now constitute a powerful constraint on theoretical galaxy and MBHBs evolution models, already excluding some of the more optimistic predictions.

In this paper, we present the analysis of MBHB populations generated in the novel semi-analytic galaxy evolution model SHARK (Lagos et al. 2018) in terms of their detectability with LISA and PTA. The model was successfully tested against standard observational data (i.e., stellar mass function, hydrogen mass function, BH–bulge mass relation, etc.) and was also used in predictions for future missions, namely, SKA and Athena (Amarantidis et al. 2019). We note that there has not been any study of the GW emission produced by MBHBs naturally occurring in the hierarchical formation scenario of SHARK and therefore we find it is important to fill in this gap. A similar research was conducted by various groups, using both semi-analytic models (SAM; Sesana et al. 2011; Dayal et al. 2019; Barausse et al. 2020) and hydrodynamical simulations (HD; Salcido et al. 2016; Kelley et al. 2017a; Katz & Larson 2019; Katz et al. 2020). Globally, the majority of these works trace a similar set of galaxy and MBHB evolution processes, however, they naturally utilize different models to describe them. Given the insufficiency of observational data and the complexity of discussed systems, we are currently unable to place any constraints that would be strong enough to rule out a substantial number of these models. Detection rates for LISA provided by previous works cover a broad range of several to several hundreds events per year (depending on the assumptions related to MBHB evolution and seeding scenarios), where generally more detections are predicted by SAMs mainly due to larger volumes (see, e.g., discussion in Katz et al. 2020; Amarantidis et al. 2019). Similarly, estimates of the GWB amplitude in the PTA range vary between 10⁻¹⁷ <  A <  10⁻¹⁴ (Rosado et al. 2015), although most works that include a modeling of MBHB evolution (e.g., due to interactions with galactic environment and orbit eccentricity; Kelley et al. 2017b; Bonetti et al. 2019; Barausse et al. 2020) predict lower amplitudes of A <  10⁻¹⁵, which are in agreement with observationally constrained upper limits. We note that a recent detection of a common spectrum process announced by all three major PTAs (CSP; Arzoumanian et al. 2020; Goncharov et al. 2021; Chen et al. 2021) raised new questions and revived ongoing discussions. Despite the fact that the amplitude of the signal is higher than previous limits (A_CSP ∼ 2 × 10⁻¹⁵), the astrophysical origin of the signal (primarily as the GWB) is neither confirmed nor ruled out. From a theoretical point of view – there is no consensus. There are studies which prove that CSP can indeed be the GWB, based on our current knowledge (Middleton et al. 2021), but on the other hand, other studies have struggled to find realistic models predicting such high GWB amplitude (Izquierdo-Villalba et al. 2022). All of the above implies a vital need for further studies.

With this work, we would like to initiate our search for more robust methods of model comparison by utilizing the open-source community code SHARK. In contrast to the majority of previous studies, our methods can be easily reproduced, modified, and updated within a well-built framework. Here, we present our first studies of SHARK capabilities within the GW regime using a few functionalities of the code and simplifying assumptions, which we will extend in a future work.

The outline of the paper is as follows. In Sect. 2, we describe SHARK and its physical prescriptions, particularly those that are important for MBH evolution. We also present the course of our calculations of LISA detectability and the GWB at nHz frequencies. In Sect. 3, we show our results based on a set of MBH populations produced by SHARK, including merger rates, mass ratios, GWB, and predictions for the LISA and PTA detection rates. We discuss the results and present our conclusions in Sects. 4 and 5, respectively. The cosmological parameters are taken from Planck Collaboration XIII (2016), specifically h = 0.6751, Ω_m = 0.3121, Ω_λ = 0.6891.

2. Methods

2.1. SHARK physics

Below, we briefly summarize all physical processes that can be currently modeled with SHARK and we refer the reader to Lagos et al. (2018) for more details. We place a particular emphasis on the prescriptions directly related to MBHs as the main focus of this paper is to study their impact on low-frequency GW emission (equations are taken directly from Lagos et al. 2018).

Before moving to semi-analytic prescriptions, we would like to describe the basis of SHARK, namely, dark matter (DM) halo catalogs and merger trees from N-body simulations. Notably, SHARK is not restricted to just one particular realization, but instead, it allows the user to evolve galaxies based on any chosen N-body DM catalog. Here, we use L210N1536 simulation from the SURFs suite (Elahi et al. 2018) with a volume of (210 cMpc h⁻¹)³ and a spatial resolution of 4.5 ckpc h⁻¹ (where cMpc and ckpc indicate comoving units). There are 1536³ simulated particles, with a single mass of 2.21 × 10⁸ M_⊙ h⁻¹. Currently, the DM halo mass resolution is 10⁹ M_⊙ and galaxy stellar mass resolution is 10⁷ M_⊙. Future releases of SURFs will aim at resolving halo and stellar masses down to 10⁶ M_⊙ and 10⁴ M_⊙, respectively, which will be of particular importance for studying poorly constrained coevolution processes in dwarf galaxies.

2.1.1. Gas cooling

Gas cooling is one of the most important processes influencing our results and we tested both of the available models, namely Croton06 (Croton et al. 2006) and Benson10 (Benson 2010). The main difference here is the way of calculating cooling timescales.

Croton06. The assumption here is that the cooling timescale t_cool is on the order of the dynamical timescale t_ff, both measured at the cooling radius (t_cool(r_cool)/t_ff(r_cool)=1). The cooling rate in the hot halo mode (when the cooling radius is smaller than the virial one) is then given by:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{cool}}=4\pi {\rho }_{\mathrm{g}}(r_{\mathrm{cool}})r_{\mathrm{cool}}^2{\dot{r}}_{\mathrm{cool}},\end{array}$$(1)

where ρ_g(r_cool) is the gas density profile.

Benson10. In this model, the cooling timescale at the cooling radius is equal to the halo age, given by:

$$\begin{array}{c}\hfill t_{\mathrm{halo}}=\frac{{\int }_0^{\mathrm{t}}[T_{\mathrm{V}}(t^{\prime })M_{\mathrm{gas}}(t^{\prime })/t_{\mathrm{cool}}(t^{\prime })]\mathrm{d}{t}^{\prime }}{T_{\mathrm{V}}(t)M_{\mathrm{gas}}(t)/t_{\mathrm{cool}}(t)},\end{array}$$(2)

where t is the current time and T_V is the virial temperature.

2.1.2. Disk instabilities

The criterion for instability was described in Ostriker & Peebles (1973) and Efstathiou et al. (1982) as:

$$\begin{array}{c}\hfill ϵ=\frac{V_{\mathrm{circ}}}{\sqrt{1.68G{M}_{\mathrm{disk}}/r_{\mathrm{disk}}}},\end{array}$$(3)

where V_circ is the maximum circular velocity, r_disk is the half-mass disc radius, and M_disk is the disc mass (gas plus stars). If ϵ <  ϵ_disk, then the disk is considered unstable and ϵ_disk is a free parameter in SHARK.

2.1.3. Galaxy mergers

Orbits of satellite galaxies decay due to dynamical friction with halo material as well as other satellites after entering central galaxy virial radius. This decay time before galaxies merge is calculated following Lacey & Cole (1993) as:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{merge}}=f_{\mathrm{df}}{\mathrm{\Theta }}_{\mathrm{orb}}{\tau }_{\mathrm{dyn}}\left(\frac{0.3722}{\ln ({\mathrm{\Lambda }}_{\mathrm{Coulomb}})}\right)\frac{M}{M_{\mathrm{sat}}},\end{array}$$(4)

where

$$\begin{array}{c}\hfill {\mathrm{\Theta }}_{\mathrm{orb}}={\left(\frac{J}{J_{\mathrm{c}}(E)}\right)}^{0.78}\left(\frac{r_{\mathrm{c}}(E)}{R_{\mathrm{V}}}\right),\end{array}$$(5)

which is a function of orbital parameters: initial angular momentum (J) and energy (E) of the satellites orbit; J_c(E) and r_c(E) correspond to angular momentum and radius but for a circular orbit with satellite energy. Then, ${\tau }_{\mathrm{dyn}}=\frac{R_{\mathrm{V}}}{V_{\mathrm{V}}}$ is the dynamical timescale of the halo (for virial radius R_V and velocity V_V), ${\mathrm{\Lambda }}_{\mathrm{Coulomb}}=ln\left(\frac{M}{M_{\mathrm{sat}}}\right)$ is the Coulomb logarithm, where M is the halo mass of the central galaxy and M_sat is the mass of the satellite, including the mass of its halo. f_df is an adjustable parameter corresponding to approximation of decay time below the resolution of the simulations (here set to 0.1).

2.1.4. BH growth and AGN feedback

Besides galaxy mergers, the growth of MBHs is assumed to occur during two phases. The first one, the so-called “quasar mode”, is characterized by accretion of material during starburst event, which can be triggered either by galaxy merger or disk instabilities. It is a rapid process resulting in high accretion rates and very efficient growth of the MBH. It can be calculated via (Kauffmann & Haehnelt 2000):

$$\begin{array}{c}\hfill \delta m_{\mathrm{BH},\mathrm{q}}=f_{\mathrm{mbh}}\frac{m_{\mathrm{gas}}}{1+{(v_{\mathrm{mbh}}/V_{\mathrm{V}})}^2},\end{array}$$(6)

where m_gas and V_V are the cold gas mass reservoir of the starburst and the virial velocity, respectively; and f_mbh = 8 × 10⁻³ is responsible for normalization of the BH–bulge relation and v_mbh = 400 km s⁻¹ regulates the binding energy of the system relative to V_V.

Accretion rate is calculated for a bulge accretion timescale (τ_bulge):

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{BH},\mathrm{q}}=\frac{\delta m_{\mathrm{BH},\mathrm{q}}}{{\tau }_{\mathrm{bulge}}}·\end{array}$$(7)

The second growth channel is accretion during so called “radio mode” and it can be modeled following Croton16 (Croton et al. 2016) or Bower06 (Bower et al. 2006). In this work, we explore the results for both of these models.

In the case of the Croton16 model, the Bondi-Hoyle like accretion mode is assumed and completed by “maximal cooling flow” model (Nulsen & Fabian 2000), giving:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{BH},\mathrm{r}}={\kappa }_{\mathrm{R}}\frac{15}{16}\pi G\mu m_{\mathrm{p}}\frac{{\kappa }_{\mathrm{B}}{T}_{\mathrm{V}}}{\mathrm{\Lambda }(T_{\mathrm{V}})}{m}_{\mathrm{BH}},\end{array}$$(8)

where κ_R is a free parameter regulating accretion rate, κ_B, Λ(T_V) are the Boltzmann constant and the cooling function depending on virial temperature, T_V, and halo metalicity, μm_p, is the mean mass per gas particle. MBH luminosity is then given by:

$$\begin{array}{c}\hfill L_{\mathrm{MBH}}=\eta {\dot{M}}_{\mathrm{BH},\mathrm{r}}{c}^2,\end{array}$$(9)

where η is the luminosity efficiency. We explore two model variants with η set to either 0.1 (10% of the infalling matter is converted into radiation) or 0.4 (maximum value for a spinning black hole, 60% of matter contributes to MBH growth and the rest is emitted as radiation).

In Bower06, cooling depends on the free-fall and cooling radii. If the free-fall radius is much smaller than the cooling radius, then gas is accreted on a free-fall timescale. Otherwise, quasi-hydrostatic cooling is assumed (it is also the case of an effective AGN feedback). For the latter, MBH growth rate is given by:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{BH},\mathrm{r}}=\frac{L_{\mathrm{cool}}}{0.2c^2},\end{array}$$(10)

where L_cool is the cooling luminosity. It is assumed that an active galactic nucleus (AGN) will quench gas cooling if the available AGN power is related to the cooling luminosity as:

$$\begin{array}{c}\hfill L_{\mathrm{cool}}<f_{\mathrm{edd}}{L}_{\mathrm{Edd}},\end{array}$$(11)

where L_Edd is the Eddington luminosity of the MBH and f_edd is a free parameter in the range of 0.0001−0.1. It allows us to manually set the threshold for effective gas quenching as a fraction of Eddington luminosity. We tested two models with f_edd equal 0.01 and 0.001, which means that in the latter case, gas cooling will be quenched at lower AGN luminosities.

It is important to note that Croton16 and Bower06 require different approaches to modeling gas cooling. The latter requires an explicit calculation and comparison of cooling and dynamical timescales and so it can be used only with the Benson10 gas cooling model, while Croton16 is used with the Croton06 model (see Sect. 2.1.1).

The major consequence of the above is that in Bower06 models, the AGN feedback will have a stronger dependence on redshift, especially in case of the most massive galaxies (when cooling radius is calculated based on the halo age, it will be shorter at earlier times). This would mean that the higher the redshift, the more probable it is that cooling will occur on free-fall timescales and in this case, AGN feedback is ineffective in Bower06 models (Bower et al. 2006). Another important difference to note is that in Cronon16, the AGN feedback depends on gas properties and MBH mass, while in Bower06 it is limited only by the Eddington luminosity (Bower et al. 2006).

2.1.5. Other

Star formation. We can choose between four different models, which either calculate molecular gas densities from galaxy properties or employ an observationally constrained value. Here, we use the model from (Blitz & Rosolowsky 2006) based on observations of nearby galaxies (Leroy et al. 2013).

Stellar and SNe feedback. Here we discriminate two components: gas outflows from galaxy M_gal and from halo M_halo. The rates of the two outflows are always related in the same way, however M_gal depends on the supernova (SN) feedback, which can be calculated with six different models. We used Lagos et al. (2013).

Reincorporation of ejected gas. The ejected gas is added back to the galaxy on a timescale that depends on their masses. This process is fully steerable, mainly by the τ_reinc scaling parameter (Henriques et al. 2013). If set to 0, an instantaneous incorporation is implemented. Here, we used τ_reinc = 25 Gyr.

Photoionization feedback. In this case, we have two models to choose from: Lacey et al. (2016) and Sobacchi & Mesinger (2013) – and we used the latter. The difference between the two is that in Sobacchi & Mesinger (2013), the condition for halo gas cooling is redshift-dependent and employs a fixed value for the UV background z_cut = 10; however, all model parameters can be adjusted freely.

Environmental effects. Currently, there are two approaches that are generally implemented to model the halo gas of the satellite galaxies. In the first case, it is assumed that the hot halo gas is stripped from the galaxy as soon as it becomes a satellite (cold gas is maintained). In the second approach, satellite galaxies retain their hot halo gas and continue to use it as long as it is available. In our models, we used the first approach.

2.2. MBH populations

We generated 12 models that vary in three main aspects: (1) halo and MBH seed masses, (2) MBH growth models, and (3) AGN feedback and BH growth parameters. The most relevant simulation parameters and their variations are summarized in Table 1.

For each simulated model, we drew up a list of galaxies that merge within the Hubble time and contain a central MBH. The primary list consists of MBH masses m1, m2, and redshift, z, of the merger, the secondary list contains detailed information on merging galaxies that will be used to future studies. The percentage of mergers of galaxies both containing a central MBH with respect to all mergers varies between ∼8% and ∼60% for models with most or least massive MBH seed and halo masses, respectively. It is also worth noting that ∼85% of all mergers in all models are minor, namely, they fulfill the condition: 0.1 <  M_s/M_p ≤ 0.3, where M_p and M_s are total masses of the primary and secondary galaxy, respectively.

2.3. Assumptions for the GW emission

The current version of SHARK does not trace the binary evolution of MBHs, which simply merge immediately after the merger of their host galaxies (we note that this time is different from the DM halo merger time because of a delay due to dynamical friction, as explained in Sect. 2.1.3). In our work we did not include any kind of post-processing (apart from a few assumptions) and we tested performance of SHARK “as is”, therefore results presented here should be considered as upper limits. Moreover, we will use these results as our control sample for further studies which will be presented elsewhere.

In order to make the first predictions for GW detectors, we have to either calculate or assume specific timescales upon which the emission is produced (it is related to the orbital frequency evolution with time – the information which we are missing from the simulations). MBHBs that could be detected with LISA will fall in the mass range of about 10³ − 10⁷ M_⊙, which corresponds to merger times (with respect to mission lifetime) from days to years. Calculation of an expected signal-to-noise ratio (S/N) requires these times to be specified, and in most recent studies, they are drawn from a uniform distribution of 0−4 years (marking the nominal length of LISA’s operation). A more realistic approach would be to also include sources merging after LISA observations but for which we could still detect the inspiral phase. We therefore use a wider distribution of merger times (see Sect. 2.3.1) limited mostly by computational expenses. We note however, that the true statistical relevance of a given merger time distribution would be revealed for models including MBHB evolution processes.

Gravitational wave background, which is a primary target of worldwide PTA efforts, is expected to be dominated by signals from the most massive MBHBs (M >  10⁸ M_⊙), for which merger times could be on the order of billions of years. In the first instance, the amplitude of the GWB would depend on galaxy merger rates and the number density of MBHBs, which is virtually unknown due to the lack of observational data. Furthermore, there is a number of processes which will further modify the spectral shape and strain of the signal that are mostly dependent on the relations between MBHs and their hosts. These include growth mechanisms, binary formation channels, and orbital separation hardening (see e.g., Kelley et al. 2017b,a). In our study we do not address these issues directly and instead assume a fixed time delay between the merger of host galaxies and their MBHs (see Sect. 2.3.2). This allows us to focus entirely on the influence of different growth and feedback models and test if current version of SHARK is in agreement with PTA observational constraints.

2.3.1. LISA

We calculated the S/N for LISA four-year long mission using gwsnrcalc Python package included in the BOWIE analysis tool (Katz & Larson 2019), which uses PhenomD for waveform generation (Husa et al. 2016; Khan et al. 2016). The advantage of the code is that it enables us to perform calculations for all phases separately (inspiral, merger, rigndown), as well as including BH spins and taking into account binaries that may merge outside LISA’s observational window. In our calculations, we follow a similar procedure as in Katz et al. (2020): we found the number of binaries merging within 50 years using Eq. (13) (without any assumptions on the S/N) and randomly chose a list of our sources from the given distribution, assigning to each of them a merger time from 0 to 50 years (with 1 year step) and start-end times of the event with respect to LISA operations. We leave a detailed spin analysis for future works and thus we set a constant value for all MBHs here, depending on the chosen luminosity efficiency η (for η = 0.1 we set a = 0.1 and η = 0.4 would correspond to maximally spinning BH with a = 1).

The averaged S/N is calculated via (Robson et al. 2019):

$$\begin{array}{c}\hfill ⟨{\rho }^2⟩=\frac{16}{5}{\displaystyle {\int }_0^{\infty }}\frac{h_{\mathrm{c}}^2}{h_{\mathrm{N}}^2}\frac{1}{f}\mathrm{d}f,\end{array}$$(12)

where h_c and h_N are the characteristic strains of the gravitational wave and detector sensitivity, respectively.

We set a detection S/N threshold to 7. The expected number of mergers per year is then calculated via (Arun et al. 2009):

$$\begin{array}{c}\hfill \frac{d^{2}N}{\mathrm{d}z\mathrm{d}t}=4\pi c{N}_{\mathrm{com}}(z){\left(\frac{d_{\mathrm{L}}(z)}{1+z}\right)}^2[{\mathrm{year}}^{-1}],\end{array}$$(13)

where d_L(z) is the luminosity distance, c is the speed of light, and

$$\begin{array}{c}\hfill N_{\mathrm{com}}(z)=\frac{d^{2}n(z)}{\mathrm{d}z\mathrm{d}{V}_{\mathrm{c}}}\approx \frac{N(z)}{\mathrm{\Delta }z{V}_{\mathrm{c}}}\end{array}$$(14)

is the BH merger rate density (per unit comoving volume per redshift interval) and V_c is the comoving volume of our simulations (Katz & Larson 2019).

2.3.2. PTA

In order to characterize the GWB from the whole population of MBHs for each of our models, we start by calculating the characteristic strain h_c following Rosado et al. (2015). Here, we limited our analysis to binaries with total masses M ≥ 10⁸ M_⊙ and z ≤ 2 because we expect only the most massive and relatively close sources to contribute to GWB at nHz frequencies available for PTA. All of these sources were then assigned with two parameters drawn from uniform distributions, namely: ι (binary inclination with respect to the line of sight) and t (estimated time to merger). In case of the latter, we probed two scenarios, with t from 0 to either 100 Myr or 1 Gyr.

The strength of GWB can be measured in terms of characteristic strain, which is an incoherent superposition of strain from each binary as a function of frequency:

$$\begin{array}{c}\hfill h_{\mathrm{c}}^2=\frac{{\sum }_k{h}_k^2{f}_k}{\mathrm{\Delta }f},\end{array}$$(15)

where summation is running over all sources at a given frequency bin determined by the whole observation time T of an array of pulsars (Δf = 1/T). The average strain for each binary is given by:

$$\begin{array}{c}\hfill h_k=A\sqrt{\frac{1}{2}(a^2+b^2)},\end{array}$$(16)

where

$$\begin{array}{c}\hfill A=2\frac{G^{5/3}{\mathcal{M}}^{5/3}{(\pi f)}^{2/3}}{c^{4}r}\end{array}$$(17)

is the dimensionless amplitude of the signal, f is the rest frame frequency, and

$$\begin{array}{cc}& a=1+{cos}^2(\iota )\hfill \\ \hfill & b=-2cos(\iota ),\hfill \end{array}$$(18)

which are two contributions from wave polarizations. The comoving distance r is defined as:

$$\begin{array}{c}\hfill r=\frac{c}{H_{100}h}{\displaystyle {\int }_0^z}{({\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^3+{\mathrm{\Omega }}_{\lambda })}^{-1/2}\mathrm{d}z.\end{array}$$(19)

Given the time to merger assigned to every binary, we calculate the frequency of its GW emission via:

$$\begin{array}{c}\hfill f={\left(\frac{8}{3}{(\pi )}^{8/3}\frac{96}{5}\frac{G^{5/3}}{c^5}{\mathcal{M}}^{5/3}t\right)}^{-3/8}.\end{array}$$(20)

3. Results

Lagos et al. (2018) presented a detailed analysis of a model, which is marked as B4H10-n1 here (BH seed mass of 10⁴ M_⊙ and halo seed mass 10¹⁰ M_⊙, Croton16 model with η = 0.1). Additionally, we extend the analysis to a wider set of MBH/halo seed masses and two different models of MBH growth as summarized in Table 1.

Figure 1 shows merger rates per year calculated with Eq. (13) for both MBH growth models (Croton16 and Bower06) and three seed scenarios. The plot does not show models with changed η and f_edd parameters as they overlap with the default realizations (they do not affect the merger rates). It can be seen that the highest merger rate for Bower06 models is shifted towards higher redshifts with respect to Croton16 models. Moreover, the shift is growing bigger with increasing seed masses, however the Croton16 and Bower06 models are indistinguishable for the lightest seed scenario (B3H9 plotted with red color). This comes from the fact that, as noted in Croton et al. (2006), Bower et al. (2006), the cooling timescale is irrelevant for small halos. Additionally, we could expect that if the two considered models use distinct relations between cooling and dynamical radii (see Sect. 2.3.1), the results will diverge for larger halos. Lastly, Croton16 models (solid) show a slightly larger difference (approximately by a factor of 1.5) in merger rates between the medium (purple) and highest (green) seed masses than Bower06 models (dashed).

Fig. 1. Merger rates per year for Croton16 and Bower06 models including all galaxies with a central MBH. The plot shows the basis MBH growth models (with η = 0.1 and fedd = 0.01) as variations in these parameters does not affect the result.

Below, we present the results of calculating low frequency GW emission from inspiraling binaries for all 12 MBH populations. In Sect. 3.1, we focus on mHz band that will be probed by LISA and in Sect. 3.2, we show the predicted GWB at nHz frequencies which is the main target of PTA observations.

3.1. LISA

Table 2 contains the list of all 12 models analysed in this work along with merger and detection rates per year for the LISA mission. In the third column, we report the percentage of mergers where both galaxies contain a central MBH with respect to all mergers. As seen in Fig. 1, the total number of mergers per year is higher for Bower06 models (except the lightest seed scenarios). This results in approximately ten more detections with respect to Croton16 models for each corresponding seed scenario.

In general, detection rates vary between ∼60 and ∼7 per year for the least and most massive seed models, respectively. In the case of the variable η parameter for Croton16 models, the ones with η = 0.4 result in a few more detections. This mostly comes from the fact that GW signal detected by LISA will be higher for systems with larger spins. In case of Bower06 models with varying f_edd, we do not find any significant differences in detection rates which means that they are purely statistical and originate from our sampling procedure.

Figure 2 shows mass ratios for all binaries (solid lines) and those detected within four years of LISA operation (dashed lines). Here, we also report only the default versions of both Croton16 and Bower06 models. From the plot we can see that both models (in all seed scenarios) favor binaries with low mass ratios (q <  10⁻²). Particularly, there is a depletion of binaries with q in the range of 10⁻² − 10⁻¹ for the least massive seed scenario (red color) however it’s less significant for Bower06 models. This results in a much smaller detection rate in this q regime. For higher mass seeds, Croton16 models show a plateau, producing a relatively even sample of binaries with various mass ratios above q = 10⁻⁴, while in case of Bower06 models, the counts are steadily growing towards higher mass ratios and slightly dropping down at the end where q = 1. As noted in Sect. 2.2, the vast majority of mergers are minor (irrespectively of the model), which explains the distribution of q in Fig. 2.

Fig. 2. Mass ratios for all binaries generated in SHARK (solid line) and those detected over four years by LISA (dashed line) for default Croton16 and Bower06 models, i.e., with η = 0.1 and fedd = 0.01.

Comparing Figs. 1 and 2 with information from Table 2, we can see that higher detection rates in case of Bower06 models may come from two main factors: (a) these models are characterized with higher merger rates (there are more halos containing a MBH); and (b) their binary population consists of systems with higher mass ratios (producing stronger signal in the detector). This would mean, that apart from being more effective in seeding halos with MBHs (due to directly comparing cooling and dynamical radii; see Sects. 2.1.1 and 2.1.4), the Bower06 model results in a more even MBH mass distribution. We find that there are approximately 19% and 45% more mergers with q >  10⁻² for Bower06 medium and high seed scenarios, respectively.

3.2. PTA

In terms of the nHz GWB, it is practically impossible to distinguish between various seed scenarios; in addition, the two BH growth models result in a comparable signal amplitude, which is shown in Fig. 3. All Bower06 models produce a slightly higher characteristic strain due to the same two factors already mentioned for LISA case (higher merger rates and mass ratios), but the difference is miniscule and thus not statistically relevant. As GWB is produced by only the most massive (M >  10⁸ M_⊙) and closest (z <  2) binaries and given that MBH evolution and growth are partially driven by mergers, the memory of the seeding scenario is lost.

Fig. 3. GWB characteristic strain comparison between default Croton16 and Bower06 models (η = 0.1 and fedd = 0.01). Panel on the left shows models with time to coalescence up to 100 Myr, while on the right, the maximum time is 1 Gyr. The black line and blue diamond show the observationally constrained limit (PPTA; Shannon et al. 2015) and estimated SKA sensitivity (at the most most sensitive frequency), respectively.

However, GWB might carry important information on MBH binary evolution after the galaxy merger. This regime is particularly intriguing as it’s poorly understood due to our limited resolution in both theory and observations. As we mentioned before, in this work, we did not focus on binary evolution and instead we qualitatively showed how the time available for binary formation and evolution affects the predicted GWB strain amplitude.

Assuming an optimistic scenario, where all galaxy mergers result in MBH binary formation and subsequently coalescence within 100 Myr we obtain a mean characteristic strain amplitude at the reference frequency (f = 1/yr) h_c, yr = 1.087 × 10⁻¹⁵ averaged over all models. When the coalescence time is additionally delayed up to 1 Gyr, the average strain amplitude drops down to h_c, yr = 1.376 × 10⁻¹⁶. Both scenarios are shown in two panels of Fig. 3 (shorter coalescence time on the left plot and longer time on the right).

What is extremely important when considering this scenario is that current PTA observations already put rigorous constraints on possible GWB – we compare them with our results directly in Table 3. Additionally, both panels of Fig. 3 show sensitivity curve of PPTA (Shannon et al. 2015) for better visibility. These results can lead to a few conclusions. First, we can see that both coalescence times we applied fit within observational constraints, however, GWB amplitude for the 100 Myr scenario is close to currently achievable sensitivities. On the other hand, a recently detected CSP has a slightly higher (on average by a factor of ∼2) amplitude which is also summarized in Table 3. Notably, this apparent discrepancy may have been caused by imperfect solar system ephemeris models used in previous studies (for a discussion see Arzoumanian et al. 2020). Nevertheless, our predicted amplitudes lie below the observed CSP, which would disfavor the intensely debated hypothesis that the signal is indeed the GWB. The second conclusion is that increasing the delay time between galaxy and MBH merger by an order or magnitude decreases the characteristic strain amplitude by approximately half an order of magnitude. Moreover, the delay time is the only parameter that significantly affects the resultant GWB in our studies. A further decrease in the amplitude may be caused by additional processes, which we have omitted here, for instance, recoil kicks and ejections (lowering the number density of MBHBs), binaries stalling at ∼kpc separations and never reaching the GW emission regime or longer delay times caused by environmental coupling. Here, we show that even considering our optimistic upper limits, current PTAs are still incapable of detecting the GWB. We note however, that given fast advances in instrumentation and analysis (along with ever increasing timing baseline), our sensitivities are increasing by approximately half an order of magnitude with every new data release (see e.g., Aggarwal et al. 2019). What is more, future detectors like SKA are expected to reach at least two orders of magnitude lower compared to today’s experiments (shown in Fig. 3 with blue diamond) and will have a significant impact on the detection and characterization of GWB.

4. Discussion

To gain more insight into our MBH population properties, in Fig. 4 we show a BH mass function (BHMF) for a default Croton16 model at four redshift values (z = 0, 1, 2, 4) and compare it with three different observational estimates. It can be seen that the mass function is reproduced reasonably well at the lowest redshift (with a mild overestimation for M_MBH >  10⁹ M_⊙) but, generally, there is a clear discrepancy between the model and observations for the least and most massive MBHs at higher redshifts. This is a common feature with regard to the majority of simulations discussed in the literature (see e.g., Habouzit et al. 2021) and the reason can be twofold. First, all of the simulations are limited by resolution and volume which directly impact their ability to reproduce the lowest and highest mass BH populations. The second issue can be noticed immediately when looking at Fig. 4. As observational BH mass estimates are mostly limited to the local Universe, deriving BHMF evolution with redshift requires adopting some assumptions or models, e.g. various MBH–galaxy scaling relations (for an extensive review see Kelly & Merloni 2012). As a result each estimated BHMF will be different, especially at highest redshifts and any comparison with simulations should be done with caution.

Fig. 4. Black hole mass function (BHMF) for z = 0, 1, 2, 4 from Croton16 model B4H10 (pink solid line). Observationally constrained mass functions are taken from Merloni & Heinz (2008) (orange circles, here z = 0.1 instead of 0), Shankar et al. (2009) (purple diamonds), and Cao (2010) (dark blue stars).

The BH mass functions for different seed scenarios and BH growth models show little difference and this could mean that the global evolution of mass assembly is similar for all of the models (shown in Fig. 5). This is due to the fact that all models are calibrated with the same M_MBH − M_bulge relation, which regulates the rate of gas inflow onto the MBH during a starburst event (quasar mode, for more, see Sect. 2.1.4), which is also the main channel of MBH growth for all SHARK models. However, the BHMF for the two models is clearly different at the lowest mass range (M <  10⁶ M_⊙) and at large redshifts (z >  6). This could mean that variances in gas cooling prescriptions affecting the second channel of MBH growth (radio mode) mainly influence the seed formation stage; namely, Bower06 predicts more MBH seeds planted into DM halos.

Fig. 5. BHMF comparison between Croton16 and Bower06 models at z = 0, 2, 6, 10.

Finally, in Fig. 6, we show total masses of binaries detected by LISA as a function of redshift (color points) overplotted on all binaries merging within 50 years (grey points). It can be clearly seen how the difference between Croton16 and Bower06 models grows with increasing MBH seed mass. In case of Bower06, mergers start occuring at earlier times then for Croton16 for all seed scenarios and the majority of detections comprises binaries with masses in the range of 10⁵ − 10⁷ M_⊙. The most distant detected mergers occur at redshifts z = 6−7 for Bower06 high seed scenario models. Notably, the mass – redshift distribution of detections is very sensitive to model variations and could be used to make preliminary evolutionary constraints.

Fig. 6. Total masses of binaries detected by LISA over the whole four-year long mission (color) and all binaries merging in 50 years (grey) as a function of redshift. Color scale shows predicted S/N of detection. Croton16 models are shown on the left side, while Bower06 on the right. MBH and halo seed masses are increasing from top to bottom.

5. Conclusions

We present an analysis of 12 MBHB populations generated with SHARK in terms of their detectability with current and future gravitational wave detectors. Specifically, we focused on merger and detection rates for the LISA mission and the predicted GWB amplitude at nHz frequencies probed by PTA experiments.

Our findings can be summarized as follows:

1. GWB amplitude mainly depends on the delay time between host galaxy and MBH merger and is rather independent of all other parameters we investigated (different seed scenarios and growth and feedback models). The GWB spectral shape will depend on the mass spectrum, orbit eccentricity or environmental effects, however, we will investigate this in more detail in future work. Here, we assumed all binary inspirals to be circular and GW driven for simplicity.

2. The predicted GWB for all discussed models and the delay time of 100 Myr is close to the detection limit when compared with the PTA observational constraints. However, both 100 Myr and 1 Gyr delay scenarios lie well below the amplitude of common-spectrum process detected by Arzoumanian et al. (2020), Goncharov et al. (2021), Chen et al. (2021). Given that our models describe an optimistic scenario, where all MBHs form bound binaries immediately after galaxy mergers and no remnants are expelled from bulges, we can conclude that our studies disfavor the hypothesis of CSP indeed being the GWB.

3. All models are in a fairly good agreement with the observationally constrained black hole mass functions (BHMF), however, we note some discrepancies for the lowest and highest mass range at redshifts of z >  1. We note, however, that such deviations are common among most simulations (due to limited resolution and volume; e.g., Habouzit et al. 2021) and also that different observational estimates vary substantially, especially for high redshifts (due to a lack of data) and thus any comparisons should be made with caution.

4. Table 4 summarizes the key features of our and previous studies. In particular, we note types of merger trees, seeding scenarios, MBHB evolution models, and LISA detection rates. Despite the substantial differences at all levels of modeling, estimated merger rates are rather consistent (e.g., due to the parameter calibration with respect to the observed scaling relations) and therefore it could be difficult to evaluate the accuracy of each model. One way to robustly estimate the contribution of each model to the detection rates would be to perform large-scale comparative studies within one consistent framework, such as SHARK. This is the undertaking we have initiated with the present work. We showed that the mass-redshift distribution of the detected sources strongly depends on the evolution models (of both MBH seeds and their growth) and given our uniform framework, this should allow us to place new constraints on galaxy evolution when compared to the actual detections we expect in the future.

5. Finally, we also find that the mass ratios of the majority of detectable binaries will be far from unity, implying the importance of developing analysis techniques for high and extreme mass ratio inspirals.

The presented results should be considered as the upper limits given the assumptions applied in our studies, namely: all galaxy mergers result in formation of a MBHB, preconceived MBH seed masses, and merger delay times. Future works will extend our analysis to the details of MBH binary evolution based on the characteristics of merging galaxies (specifically related to the problem of merger time delays) and the possible electromagnetic signatures of MBH binary inspirals.

Acknowledgments

T.B. was suported by the Foundation for Polish Science grant TEAM/2016-3/19. We are grateful to the SHARK team for helpful discussions.

References

All Tables

All Figures

	Fig. 1. Merger rates per year for Croton16 and Bower06 models including all galaxies with a central MBH. The plot shows the basis MBH growth models (with η = 0.1 and fedd = 0.01) as variations in these parameters does not affect the result.
In the text

	Fig. 2. Mass ratios for all binaries generated in SHARK (solid line) and those detected over four years by LISA (dashed line) for default Croton16 and Bower06 models, i.e., with η = 0.1 and fedd = 0.01.
In the text

Fig. 3. GWB characteristic strain comparison between default Croton16 and Bower06 models (η = 0.1 and fedd = 0.01). Panel on the left shows models with time to coalescence up to 100 Myr, while on the right, the maximum time is 1 Gyr. The black line and blue diamond show the observationally constrained limit (PPTA; Shannon et al. 2015) and estimated SKA sensitivity (at the most most sensitive frequency), respectively.

In the text

	Fig. 4. Black hole mass function (BHMF) for z = 0, 1, 2, 4 from Croton16 model B4H10 (pink solid line). Observationally constrained mass functions are taken from Merloni & Heinz (2008) (orange circles, here z = 0.1 instead of 0), Shankar et al. (2009) (purple diamonds), and Cao (2010) (dark blue stars).
In the text

	Fig. 5. BHMF comparison between Croton16 and Bower06 models at z = 0, 2, 6, 10.
In the text

	Fig. 6. Total masses of binaries detected by LISA over the whole four-year long mission (color) and all binaries merging in 50 years (grey) as a function of redshift. Color scale shows predicted S/N of detection. Croton16 models are shown on the left side, while Bower06 on the right. MBH and halo seed masses are increasing from top to bottom.
In the text